home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / msdos / fractal / frasr182 / fractals.c < prev    next >
C/C++ Source or Header  |  1993-04-16  |  105KB  |  3,870 lines
  1. /*
  2. FRACTALS.C, FRACTALP.C and CALCFRAC.C actually calculate the fractal
  3. images (well, SOMEBODY had to do it!).  The modules are set up so that
  4. all logic that is independent of any fractal-specific code is in
  5. CALCFRAC.C, the code that IS fractal-specific is in FRACTALS.C, and the
  6. struscture that ties (we hope!) everything together is in FRACTALP.C.
  7. Original author Tim Wegner, but just about ALL the authors have
  8. contributed SOME code to this routine at one time or another, or
  9. contributed to one of the many massive restructurings.
  10.  
  11. The Fractal-specific routines are divided into three categories:
  12.  
  13. 1. Routines that are called once-per-orbit to calculate the orbit
  14.    value. These have names like "XxxxFractal", and their function
  15.    pointers are stored in fractalspecific[fractype].orbitcalc. EVERY
  16.    new fractal type needs one of these. Return 0 to continue iterations,
  17.    1 if we're done. Results for integer fractals are left in 'lnew.x' and
  18.    'lnew.y', for floating point fractals in 'new.x' and 'new.y'.
  19.  
  20. 2. Routines that are called once per pixel to set various variables
  21.    prior to the orbit calculation. These have names like xxx_per_pixel
  22.    and are fairly generic - chances are one is right for your new type.
  23.    They are stored in fractalspecific[fractype].per_pixel.
  24.  
  25. 3. Routines that are called once per screen to set various variables.
  26.    These have names like XxxxSetup, and are stored in
  27.    fractalspecific[fractype].per_image.
  28.  
  29. 4. The main fractal routine. Usually this will be StandardFractal(),
  30.    but if you have written a stand-alone fractal routine independent
  31.    of the StandardFractal mechanisms, your routine name goes here,
  32.    stored in fractalspecific[fractype].calctype.per_image.
  33.  
  34. Adding a new fractal type should be simply a matter of adding an item
  35. to the 'fractalspecific' structure, writing (or re-using one of the existing)
  36. an appropriate setup, per_image, per_pixel, and orbit routines.
  37.  
  38. --------------------------------------------------------------------   */
  39.  
  40.  
  41. #include <stdio.h>
  42. #include <stdlib.h>
  43. #include <float.h>
  44. #include <limits.h>
  45. #include <string.h>
  46. #ifdef __TURBOC__
  47. #include <alloc.h>
  48. #elif !defined(__386BSD__)
  49. #include <malloc.h>
  50. #endif
  51. #include "fractint.h"
  52. #include "mpmath.h"
  53. #include "helpdefs.h"
  54. #include "fractype.h"
  55. #include "prototyp.h"
  56.  
  57. #define NEWTONDEGREELIMIT  100
  58.  
  59. extern int using_jiim;
  60. extern _CMPLX  initorbit;
  61. extern _LCMPLX linitorbit;
  62. extern char useinitorbit;
  63. extern double fgLimit;
  64. extern int distest;
  65.  
  66. #define CMPLXsqr_old(out)    \
  67.    (out).y = (old.x+old.x) * old.y;\
  68.    (out).x = tempsqrx - tempsqry
  69.  
  70. #define CMPLXpwr(arg1,arg2,out)   (out)= ComplexPower((arg1), (arg2))
  71. #define CMPLXmult1(arg1,arg2,out)    Arg2->d = (arg1); Arg1->d = (arg2);\
  72.      dStkMul(); Arg1++; Arg2++; (out) = Arg2->d
  73.  
  74. #define CMPLXmult(arg1,arg2,out)  \
  75.     {\
  76.        _CMPLX TmP;\
  77.        TmP.x = (arg1).x*(arg2).x - (arg1).y*(arg2).y;\
  78.        TmP.y = (arg1).x*(arg2).y + (arg1).y*(arg2).x;\
  79.        (out) = TmP;\
  80.      }
  81. #if 0
  82. #define CMPLXmult(arg1,arg2,out)  {CMPLXmult2(arg1,arg2,&out);}
  83. #endif
  84.  
  85. void CMPLXmult2(_CMPLX arg1,_CMPLX arg2, _CMPLX *out)
  86. {
  87.     _CMPLX TmP;
  88.     TmP.x = (arg1).x*(arg2).x - (arg1).y*(arg2).y;
  89.     TmP.y = (arg1).x*(arg2).y + (arg1).y*(arg2).x;
  90.     *out = TmP;
  91. }
  92.  
  93. #define CMPLXadd(arg1,arg2,out)    \
  94.     (out).x = (arg1).x + (arg2).x; (out).y = (arg1).y + (arg2).y
  95. #define CMPLXsub(arg1,arg2,out)    \
  96.     (out).x = (arg1).x - (arg2).x; (out).y = (arg1).y - (arg2).y
  97. #define CMPLXtimesreal(arg,real,out)   \
  98.     (out).x = (arg).x*(real);\
  99.     (out).y = (arg).y*(real)
  100.  
  101. #define CMPLXrecip(arg,out)    \
  102.    { double denom; denom = sqr((arg).x) + sqr((arg).y);\
  103.      if(denom==0.0) {(out).x = 1.0e10;(out).y = 1.0e10;}else\
  104.     { (out).x =  (arg).x/denom;\
  105.      (out).y = -(arg).y/denom;}}
  106.  
  107. extern int xshift, yshift;
  108. extern int biomorph;
  109. extern int forcesymmetry;
  110. extern int symmetry;
  111. _LCMPLX lcoefficient,lold,lnew,lparm, linit,ltmp,ltmp2,lparm2;
  112. long ltempsqrx,ltempsqry;
  113. extern int decomp[];
  114. extern double param[];
  115. extern int potflag;                   /* potential enabled? */
  116. extern double f_radius,f_xcenter,f_ycenter;    /* inversion radius, center */
  117. extern double xxmin,xxmax,yymin,yymax;           /* corners */
  118. extern int overflow;
  119. extern int integerfractal;    /* TRUE if fractal uses integer math */
  120.  
  121. int maxcolor;
  122. int root, degree,basin;
  123. double floatmin,floatmax;
  124. double roverd, d1overd, threshold;
  125. _CMPLX tmp2;
  126. extern _CMPLX init,tmp,old,new,saved,jbcfp;
  127. _CMPLX coefficient;
  128. _CMPLX    staticroots[16]; /* roots array for degree 16 or less */
  129. _CMPLX    *roots = staticroots;
  130. struct MPC    *MPCroots;
  131. extern int color, row, col;
  132. extern int invert;
  133. extern double far *dx0, far *dy0;
  134. extern double far *dx1, far *dy1;
  135. long FgHalf;
  136.  
  137. _CMPLX one;
  138. _CMPLX pwr;
  139. _CMPLX Coefficient;
  140.  
  141. extern int    colors;             /* maximum colors available */
  142. extern int    inside;             /* "inside" color to use    */
  143. extern int    outside;            /* "outside" color to use   */
  144. extern int    finattract;
  145. extern int    fractype;            /* fractal type */
  146. extern int    debugflag;            /* for debugging purposes */
  147.  
  148. extern double    param[];        /* parameters */
  149. extern long    far *lx0, far *ly0;        /* X, Y points */
  150. extern long    far *lx1, far *ly1;        /* X, Y points */
  151. extern long    delx,dely;            /* X, Y increments */
  152. extern long    delmin;             /* min(max(dx),max(dy) */
  153. extern double    ddelmin;            /* min(max(dx),max(dy) */
  154. extern long    fudge;                /* fudge factor (2**n) */
  155. extern int    bitshift;            /* bit shift for fudge */
  156. int    bitshiftless1;            /* bit shift less 1 */
  157.  
  158. #ifndef sqr
  159. #define sqr(x) ((x)*(x))
  160. #endif
  161.  
  162. #ifndef lsqr
  163. #define lsqr(x) (multiply((x),(x),bitshift))
  164. #endif
  165.  
  166. #define modulus(z)     (sqr((z).x)+sqr((z).y))
  167. #define conjugate(pz)    ((pz)->y = 0.0 - (pz)->y)
  168. #define distance(z1,z2)  (sqr((z1).x-(z2).x)+sqr((z1).y-(z2).y))
  169. #define pMPsqr(z) (*pMPmul((z),(z)))
  170. #define MPdistance(z1,z2)  (*pMPadd(pMPsqr(*pMPsub((z1).x,(z2).x)),pMPsqr(*pMPsub((z1).y,(z2).y))))
  171.  
  172. double twopi = PI*2.0;
  173. int c_exp;
  174.  
  175.  
  176. /* These are local but I don't want to pass them as parameters */
  177. _CMPLX lambda;
  178. extern double magnitude, rqlim, rqlim2;
  179. _CMPLX parm,parm2;
  180. _CMPLX *floatparm;
  181. _LCMPLX *longparm; /* used here and in jb.c */
  182. extern int (*calctype)();
  183. extern unsigned long lm;        /* magnitude limit (CALCMAND) */
  184.  
  185. /* -------------------------------------------------------------------- */
  186. /*        These variables are external for speed's sake only      */
  187. /* -------------------------------------------------------------------- */
  188.  
  189. double sinx,cosx,sinhx,coshx;
  190. double siny,cosy,sinhy,coshy;
  191. double tmpexp;
  192. double tempsqrx,tempsqry,tempsqrz,tempsqrt;
  193.  
  194. double foldxinitx,foldyinity,foldxinity,foldyinitx;
  195. long oldxinitx,oldyinity,oldxinity,oldyinitx;
  196. long longtmp;
  197. extern long lmagnitud, llimit, llimit2, l16triglim;
  198. extern periodicitycheck;
  199. extern char floatflag;
  200.  
  201. /* These are for quaternions */
  202. double qc,qci,qcj,qck;
  203.  
  204. /* temporary variables for trig use */
  205. long lcosx, lcoshx, lsinx, lsinhx;
  206. long lcosy, lcoshy, lsiny, lsinhy;
  207.  
  208. /*
  209. **  details of finite attractors (required for Magnet Fractals)
  210. **  (can also be used in "coloring in" the lakes of Julia types)
  211. */
  212. extern          int      attractors; /* number of finite attractors   */
  213. extern _CMPLX  attr[];       /* finite attractor values (f.p) */
  214. extern _LCMPLX lattr[];      /* finite attractor values (int) */
  215. extern int    attrperiod[]; /* finite attractor period */
  216.  
  217. /*
  218. **  pre-calculated values for fractal types Magnet2M & Magnet2J
  219. */
  220. _CMPLX    T_Cm1;          /* 3 * (floatparm - 1)            */
  221. _CMPLX    T_Cm2;          /* 3 * (floatparm - 2)            */
  222. _CMPLX    T_Cm1Cm2;     /* (floatparm - 1) * (floatparm - 2) */
  223.  
  224. void FloatPreCalcMagnet2() /* precalculation for Magnet2 (M & J) for speed */
  225.   {
  226.     T_Cm1.x = floatparm->x - 1.0;   T_Cm1.y = floatparm->y;
  227.     T_Cm2.x = floatparm->x - 2.0;   T_Cm2.y = floatparm->y;
  228.     T_Cm1Cm2.x = (T_Cm1.x * T_Cm2.x) - (T_Cm1.y * T_Cm2.y);
  229.     T_Cm1Cm2.y = (T_Cm1.x * T_Cm2.y) + (T_Cm1.y * T_Cm2.x);
  230.     T_Cm1.x += T_Cm1.x + T_Cm1.x;   T_Cm1.y += T_Cm1.y + T_Cm1.y;
  231.     T_Cm2.x += T_Cm2.x + T_Cm2.x;   T_Cm2.y += T_Cm2.y + T_Cm2.y;
  232.   }
  233.  
  234.  
  235. /* -------------------------------------------------------------------- */
  236. /*        Bailout Routines Macros                                                 */
  237. /* -------------------------------------------------------------------- */
  238.  
  239. static int near floatbailout()
  240. {
  241.    if ( ( magnitude = ( tempsqrx=sqr(new.x) )
  242.             + ( tempsqry=sqr(new.y) ) ) >= rqlim ) return(1);
  243.    old = new;
  244.    return(0);
  245. }
  246.  
  247. /* longbailout() is equivalent to next */
  248. #define LONGBAILOUT()    \
  249.    ltempsqrx = lsqr(lnew.x); ltempsqry = lsqr(lnew.y);\
  250.    lmagnitud = ltempsqrx + ltempsqry;\
  251.    if (lmagnitud >= llimit || lmagnitud < 0 || labs(lnew.x) > llimit2\
  252.      || labs(lnew.y) > llimit2 || overflow) \
  253.            { overflow=0;return(1);}\
  254.    lold = lnew;
  255.  
  256. #define FLOATTRIGBAILOUT()  \
  257.    if (fabs(old.y) >= rqlim2) return(1);
  258.  
  259. #define LONGTRIGBAILOUT()  \
  260.    if(labs(lold.y) >= llimit2 || overflow) { overflow=0;return(1);}
  261.  
  262. #define LONGXYTRIGBAILOUT()  \
  263.    if(labs(lold.x) >= llimit2 || labs(lold.y) >= llimit2 || overflow)\
  264.     { overflow=0;return(1);}
  265.  
  266. #define FLOATXYTRIGBAILOUT()  \
  267.    if (fabs(old.x) >= rqlim2 || fabs(old.y) >= rqlim2) return(1);
  268.  
  269. #define FLOATHTRIGBAILOUT()  \
  270.    if (fabs(old.x) >= rqlim2) return(1);
  271.  
  272. #define LONGHTRIGBAILOUT()  \
  273.    if(labs(lold.x) >= llimit2 || overflow) { overflow=0;return(1);}
  274.  
  275. #define TRIG16CHECK(X)    \
  276.       if(labs((X)) > l16triglim || overflow) { overflow=0;return(1);}
  277.  
  278. #define FLOATEXPBAILOUT()  \
  279.    if (fabs(old.y) >= 1.0e8) return(1);\
  280.    if (fabs(old.x) >= 6.4e2) return(1);
  281.  
  282. #define LONGEXPBAILOUT()  \
  283.    if (labs(lold.y) >= (1000L<<bitshift)) return(1);\
  284.    if (labs(lold.x) >=      (8L<<bitshift)) return(1);
  285.  
  286. #if 0
  287. /* this define uses usual trig instead of fast trig */
  288. #define FPUsincos(px,psinx,pcosx) \
  289.    *(psinx) = sin(*(px));\
  290.    *(pcosx) = cos(*(px));
  291.  
  292. #define FPUsinhcosh(px,psinhx,pcoshx) \
  293.    *(psinhx) = sinh(*(px));\
  294.    *(pcoshx) = cosh(*(px));
  295. #endif
  296.  
  297. #define LTRIGARG(X)    \
  298.    if(labs((X)) > l16triglim)\
  299.    {\
  300.       double tmp;\
  301.       tmp = (X);\
  302.       tmp /= fudge;\
  303.       tmp = fmod(tmp,twopi);\
  304.       tmp *= fudge;\
  305.       (X) = tmp;\
  306.    }\
  307.  
  308. static int near Halleybailout()
  309. {
  310.    if ( fabs(modulus(new)-modulus(old)) < parm2.x)
  311.       return(1);
  312.    old = new;
  313.    return(0);
  314. }
  315.  
  316. #ifndef XFRACT
  317. #define MPCmod(m) (*pMPadd(*pMPmul((m).x, (m).x), *pMPmul((m).y, (m).y)))
  318. struct MPC mpcold, mpcnew, mpctmp, mpctmp1;
  319. struct MP mptmpparm2x;
  320.  
  321. static int near MPCHalleybailout()
  322. {
  323. static struct MP mptmpbailout;
  324.    mptmpbailout = *MPabs(*pMPsub(MPCmod(mpcnew), MPCmod(mpcold)));
  325.    if (pMPcmp(mptmpbailout, mptmpparm2x) < 0)
  326.       return(1);
  327.    mpcold = mpcnew;
  328.    return(0);
  329. }
  330. #endif
  331.  
  332. /* -------------------------------------------------------------------- */
  333. /*        Fractal (once per iteration) routines            */
  334. /* -------------------------------------------------------------------- */
  335. static double xt, yt, t2;
  336.  
  337. /* Raise complex number (base) to the (exp) power, storing the result
  338. ** in complex (result).
  339. */
  340. void cpower(_CMPLX *base, int exp, _CMPLX *result)
  341. {
  342.     if (exp<0) {
  343.     cpower(base,-exp,result);
  344.     CMPLXrecip(*result,*result);
  345.     return;
  346.     }
  347.  
  348.     xt = base->x;   yt = base->y;
  349.  
  350.     if (exp & 1)
  351.     {
  352.        result->x = xt;
  353.        result->y = yt;
  354.     }
  355.     else
  356.     {
  357.        result->x = 1.0;
  358.        result->y = 0.0;
  359.     }
  360.  
  361.     exp >>= 1;
  362.     while (exp)
  363.     {
  364.     t2 = xt * xt - yt * yt;
  365.     yt = 2 * xt * yt;
  366.     xt = t2;
  367.  
  368.     if (exp & 1)
  369.     {
  370.         t2 = xt * result->x - yt * result->y;
  371.         result->y = result->y * xt + yt * result->x;
  372.         result->x = t2;
  373.     }
  374.     exp >>= 1;
  375.     }
  376. }
  377. /* long version */
  378. static long lxt, lyt, lt2;
  379.  
  380. lcpower(_LCMPLX *base, int exp, _LCMPLX *result, int bitshift)
  381. {
  382.     static long maxarg;
  383.     maxarg = 64L<<bitshift;
  384.  
  385.     overflow = 0;
  386.     lxt = base->x;   lyt = base->y;
  387.  
  388.     if (exp & 1)
  389.     {
  390.        result->x = lxt;
  391.        result->y = lyt;
  392.     }
  393.     else
  394.     {
  395.        result->x = 1L<<bitshift;
  396.        result->y = 0L;
  397.     }
  398.  
  399.     exp >>= 1;
  400.     while (exp)
  401.     {
  402.     /*
  403.     if(labs(lxt) >= maxarg || labs(lyt) >= maxarg)
  404.        return(-1);
  405.     */
  406.     lt2 = multiply(lxt, lxt, bitshift) - multiply(lyt,lyt,bitshift);
  407.     lyt = multiply(lxt,lyt,bitshiftless1);
  408.     if(overflow)
  409.        return(overflow);
  410.     lxt = lt2;
  411.  
  412.     if (exp & 1)
  413.     {
  414.         lt2 = multiply(lxt,result->x, bitshift) - multiply(lyt,result->y,bitshift);
  415.         result->y = multiply(result->y,lxt,bitshift) + multiply(lyt,result->x,bitshift);
  416.         result->x = lt2;
  417.     }
  418.     exp >>= 1;
  419.     }
  420.     if(result->x == 0 && result->y == 0)
  421.        overflow = 1;
  422.     return(overflow);
  423. }
  424. #if 0
  425. z_to_the_z(_CMPLX *z, _CMPLX *out)
  426. {
  427.     static _CMPLX tmp1,tmp2;
  428.     /* raises complex z to the z power */
  429.     int errno_xxx;
  430.     errno_xxx = 0;
  431.  
  432.     if(fabs(z->x) < DBL_EPSILON) return(-1);
  433.  
  434.     /* log(x + iy) = 1/2(log(x*x + y*y) + i(arc_tan(y/x)) */
  435.     tmp1.x = .5*log(sqr(z->x)+sqr(z->y));
  436.  
  437.     /* the fabs in next line added to prevent discontinuity in image */
  438.     tmp1.y = atan(fabs(z->y/z->x));
  439.  
  440.     /* log(z)*z */
  441.     tmp2.x = tmp1.x * z->x - tmp1.y * z->y;
  442.     tmp2.y = tmp1.x * z->y + tmp1.y * z->x;
  443.  
  444.     /* z*z = e**(log(z)*z) */
  445.     /* e**(x + iy) =  e**x * (cos(y) + isin(y)) */
  446.  
  447.     tmpexp = exp(tmp2.x);
  448.  
  449.     FPUsincos(&tmp2.y,&siny,&cosy);
  450.     out->x = tmpexp*cosy;
  451.     out->y = tmpexp*siny;
  452.     return(errno_xxx);
  453. }
  454. #endif
  455.  
  456. int complex_div(_CMPLX arg1,_CMPLX arg2,_CMPLX *pz);
  457. int complex_mult(_CMPLX arg1,_CMPLX arg2,_CMPLX *pz);
  458.  
  459. #ifdef XFRACT /* fractint uses the NewtonFractal2 code in newton.asm */
  460.  
  461. /* Distance of complex z from unit circle */
  462. #define DIST1(z) (((z).x-1.0)*((z).x-1.0)+((z).y)*((z).y))
  463. #define LDIST1(z) (lsqr((((z).x)-fudge)) + lsqr(((z).y)))
  464.  
  465.  
  466. int NewtonFractal2()
  467. {
  468.     static char start=1;
  469.     if(start)
  470.     {
  471.        start = 0;
  472.     }
  473.     cpower(&old, degree-1, &tmp);
  474.     complex_mult(tmp, old, &new);
  475.  
  476.     if (DIST1(new) < threshold)
  477.     {
  478.        if(fractype==NEWTBASIN || fractype==MPNEWTBASIN)
  479.        {
  480.       int tmpcolor;
  481.       int i;
  482.       tmpcolor = -1;
  483.       /* this code determines which degree-th root of root the
  484.          Newton formula converges to. The roots of a 1 are
  485.          distributed on a circle of radius 1 about the origin. */
  486.       for(i=0;i<degree;i++)
  487.          /* color in alternating shades with iteration according to
  488.         which root of 1 it converged to */
  489.           if(distance(roots[i],old) < threshold)
  490.           {
  491.           if (basin==2) {
  492.                   tmpcolor = 1+(i&7)+((color&1)<<3);
  493.           } else {
  494.               tmpcolor = 1+i;
  495.           }
  496.           break;
  497.           }
  498.        if(tmpcolor == -1)
  499.           color = maxcolor;
  500.        else
  501.           color = tmpcolor;
  502.        }
  503.        return(1);
  504.     }
  505.     new.x = d1overd * new.x + roverd;
  506.     new.y *= d1overd;
  507.  
  508.     /* Watch for divide underflow */
  509.     if ((t2 = tmp.x * tmp.x + tmp.y * tmp.y) < FLT_MIN)
  510.       return(1);
  511.     else
  512.     {
  513.     t2 = 1.0 / t2;
  514.     old.x = t2 * (new.x * tmp.x + new.y * tmp.y);
  515.     old.y = t2 * (new.y * tmp.x - new.x * tmp.y);
  516.     }
  517.     return(0);
  518. }
  519.  
  520.  
  521.  
  522. #endif /* newton code only used by xfractint */
  523.  
  524. complex_mult(arg1,arg2,pz)
  525. _CMPLX arg1,arg2,*pz;
  526. {
  527.    pz->x = arg1.x*arg2.x - arg1.y*arg2.y;
  528.    pz->y = arg1.x*arg2.y+arg1.y*arg2.x;
  529.    return(0);
  530. }
  531.  
  532.  
  533. complex_div(numerator,denominator,pout)
  534. _CMPLX numerator,denominator,*pout;
  535. {
  536.    double mod;
  537.    if((mod = modulus(denominator)) < FLT_MIN)
  538.       return(1);
  539.    conjugate(&denominator);
  540.    complex_mult(numerator,denominator,pout);
  541.    pout->x = pout->x/mod;
  542.    pout->y = pout->y/mod;
  543.    return(0);
  544. }
  545.  
  546. #ifndef XFRACT
  547. struct MP mproverd, mpd1overd, mpthreshold,sqrmpthreshold;
  548. struct MP mpt2;
  549. struct MP mpone;
  550. extern struct MPC MPCone;
  551. extern int MPOverflow;
  552. #endif
  553.  
  554. int MPCNewtonFractal()
  555. {
  556. #ifndef XFRACT
  557.     MPOverflow = 0;
  558.     mpctmp   = MPCpow(mpcold,degree-1);
  559.  
  560.     mpcnew.x = *pMPsub(*pMPmul(mpctmp.x,mpcold.x),*pMPmul(mpctmp.y,mpcold.y));
  561.     mpcnew.y = *pMPadd(*pMPmul(mpctmp.x,mpcold.y),*pMPmul(mpctmp.y,mpcold.x));
  562.     mpctmp1.x = *pMPsub(mpcnew.x, MPCone.x);
  563.     mpctmp1.y = *pMPsub(mpcnew.y, MPCone.y);
  564.     if(pMPcmp(MPCmod(mpctmp1),mpthreshold)< 0)
  565.     {
  566.       if(fractype==MPNEWTBASIN)
  567.       {
  568.      int tmpcolor;
  569.      int i;
  570.      tmpcolor = -1;
  571.      for(i=0;i<degree;i++)
  572.          if(pMPcmp(MPdistance(MPCroots[i],mpcold),mpthreshold) < 0)
  573.          {
  574.         if(basin==2)
  575.            tmpcolor = 1+(i&7) + ((color&1)<<3);
  576.         else
  577.            tmpcolor = 1+i;
  578.             break;
  579.          }
  580.       if(tmpcolor == -1)
  581.          color = maxcolor;
  582.       else
  583.          color = tmpcolor;
  584.       }
  585.        return(1);
  586.     }
  587.  
  588.     mpcnew.x = *pMPadd(*pMPmul(mpd1overd,mpcnew.x),mproverd);
  589.     mpcnew.y = *pMPmul(mpcnew.y,mpd1overd);
  590.     mpt2 = MPCmod(mpctmp);
  591.     mpt2 = *pMPdiv(mpone,mpt2);
  592.     mpcold.x = *pMPmul(mpt2,(*pMPadd(*pMPmul(mpcnew.x,mpctmp.x),*pMPmul(mpcnew.y,mpctmp.y))));
  593.     mpcold.y = *pMPmul(mpt2,(*pMPsub(*pMPmul(mpcnew.y,mpctmp.x),*pMPmul(mpcnew.x,mpctmp.y))));
  594.     new.x = *pMP2d(mpcold.x);
  595.     new.y = *pMP2d(mpcold.y);
  596.     return(MPOverflow);
  597. #endif
  598. }
  599.  
  600.  
  601. Barnsley1Fractal()
  602. {
  603. #ifndef XFRACT
  604.    /* Barnsley's Mandelbrot type M1 from "Fractals
  605.    Everywhere" by Michael Barnsley, p. 322 */
  606.  
  607.    /* calculate intermediate products */
  608.    oldxinitx   = multiply(lold.x, longparm->x, bitshift);
  609.    oldyinity   = multiply(lold.y, longparm->y, bitshift);
  610.    oldxinity   = multiply(lold.x, longparm->y, bitshift);
  611.    oldyinitx   = multiply(lold.y, longparm->x, bitshift);
  612.    /* orbit calculation */
  613.    if(lold.x >= 0)
  614.    {
  615.       lnew.x = (oldxinitx - longparm->x - oldyinity);
  616.       lnew.y = (oldyinitx - longparm->y + oldxinity);
  617.    }
  618.    else
  619.    {
  620.       lnew.x = (oldxinitx + longparm->x - oldyinity);
  621.       lnew.y = (oldyinitx + longparm->y + oldxinity);
  622.    }
  623.    return(longbailout());
  624. #endif
  625. }
  626.  
  627. Barnsley1FPFractal()
  628. {
  629.    /* Barnsley's Mandelbrot type M1 from "Fractals
  630.    Everywhere" by Michael Barnsley, p. 322 */
  631.    /* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
  632.  
  633.    /* calculate intermediate products */
  634.    foldxinitx = old.x * floatparm->x;
  635.    foldyinity = old.y * floatparm->y;
  636.    foldxinity = old.x * floatparm->y;
  637.    foldyinitx = old.y * floatparm->x;
  638.    /* orbit calculation */
  639.    if(old.x >= 0)
  640.    {
  641.       new.x = (foldxinitx - floatparm->x - foldyinity);
  642.       new.y = (foldyinitx - floatparm->y + foldxinity);
  643.    }
  644.    else
  645.    {
  646.       new.x = (foldxinitx + floatparm->x - foldyinity);
  647.       new.y = (foldyinitx + floatparm->y + foldxinity);
  648.    }
  649.    return(floatbailout());
  650. }
  651.  
  652. Barnsley2Fractal()
  653. {
  654. #ifndef XFRACT
  655.    /* An unnamed Mandelbrot/Julia function from "Fractals
  656.    Everywhere" by Michael Barnsley, p. 331, example 4.2 */
  657.    /* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
  658.  
  659.    /* calculate intermediate products */
  660.    oldxinitx   = multiply(lold.x, longparm->x, bitshift);
  661.    oldyinity   = multiply(lold.y, longparm->y, bitshift);
  662.    oldxinity   = multiply(lold.x, longparm->y, bitshift);
  663.    oldyinitx   = multiply(lold.y, longparm->x, bitshift);
  664.  
  665.    /* orbit calculation */
  666.    if(oldxinity + oldyinitx >= 0)
  667.    {
  668.       lnew.x = oldxinitx - longparm->x - oldyinity;
  669.       lnew.y = oldyinitx - longparm->y + oldxinity;
  670.    }
  671.    else
  672.    {
  673.       lnew.x = oldxinitx + longparm->x - oldyinity;
  674.       lnew.y = oldyinitx + longparm->y + oldxinity;
  675.    }
  676.    return(longbailout());
  677. #endif
  678. }
  679.  
  680. Barnsley2FPFractal()
  681. {
  682.    /* An unnamed Mandelbrot/Julia function from "Fractals
  683.    Everywhere" by Michael Barnsley, p. 331, example 4.2 */
  684.  
  685.    /* calculate intermediate products */
  686.    foldxinitx = old.x * floatparm->x;
  687.    foldyinity = old.y * floatparm->y;
  688.    foldxinity = old.x * floatparm->y;
  689.    foldyinitx = old.y * floatparm->x;
  690.  
  691.    /* orbit calculation */
  692.    if(foldxinity + foldyinitx >= 0)
  693.    {
  694.       new.x = foldxinitx - floatparm->x - foldyinity;
  695.       new.y = foldyinitx - floatparm->y + foldxinity;
  696.    }
  697.    else
  698.    {
  699.       new.x = foldxinitx + floatparm->x - foldyinity;
  700.       new.y = foldyinitx + floatparm->y + foldxinity;
  701.    }
  702.    return(floatbailout());
  703. }
  704.  
  705. JuliaFractal()
  706. {
  707. #ifndef XFRACT
  708.    /* used for C prototype of fast integer math routines for classic
  709.       Mandelbrot and Julia */
  710.    lnew.x  = ltempsqrx - ltempsqry + longparm->x;
  711.    lnew.y = multiply(lold.x, lold.y, bitshiftless1) + longparm->y;
  712.    return(longbailout());
  713. #elif !defined(__386BSD__)
  714.    fprintf(stderr,"JuliaFractal called\n");
  715.    exit(-1);
  716. #endif
  717. }
  718.  
  719. JuliafpFractal()
  720. {
  721.    /* floating point version of classical Mandelbrot/Julia */
  722.    /* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
  723.    new.x = tempsqrx - tempsqry + floatparm->x;
  724.    new.y = 2.0 * old.x * old.y + floatparm->y;
  725.    return(floatbailout());
  726. }
  727.  
  728. static double f(double x, double y)
  729. {
  730.    return(x - x*y);
  731. }
  732.  
  733. static double g(double x, double y)
  734. {
  735.    return(-y + x*y);
  736. }
  737.  
  738. LambdaFPFractal()
  739. {
  740.    /* variation of classical Mandelbrot/Julia */
  741.    /* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
  742.  
  743.    tempsqrx = old.x - tempsqrx + tempsqry;
  744.    tempsqry = -(old.y * old.x);
  745.    tempsqry += tempsqry + old.y;
  746.  
  747.    new.x = floatparm->x * tempsqrx - floatparm->y * tempsqry;
  748.    new.y = floatparm->x * tempsqry + floatparm->y * tempsqrx;
  749.    return(floatbailout());
  750. }
  751.  
  752. LambdaFractal()
  753. {
  754. #ifndef XFRACT
  755.    /* variation of classical Mandelbrot/Julia */
  756.  
  757.    /* in complex math) temp = Z * (1-Z) */
  758.    ltempsqrx = lold.x - ltempsqrx + ltempsqry;
  759.    ltempsqry = lold.y
  760.          - multiply(lold.y, lold.x, bitshiftless1);
  761.    /* (in complex math) Z = Lambda * Z */
  762.    lnew.x = multiply(longparm->x, ltempsqrx, bitshift)
  763.     - multiply(longparm->y, ltempsqry, bitshift);
  764.    lnew.y = multiply(longparm->x, ltempsqry, bitshift)
  765.     + multiply(longparm->y, ltempsqrx, bitshift);
  766.    return(longbailout());
  767. #endif
  768. }
  769.  
  770. SierpinskiFractal()
  771. {
  772. #ifndef XFRACT
  773.    /* following code translated from basic - see "Fractals
  774.    Everywhere" by Michael Barnsley, p. 251, Program 7.1.1 */
  775.    lnew.x = (lold.x << 1);        /* new.x = 2 * old.x  */
  776.    lnew.y = (lold.y << 1);        /* new.y = 2 * old.y  */
  777.    if(lold.y > ltmp.y)    /* if old.y > .5 */
  778.       lnew.y = lnew.y - ltmp.x; /* new.y = 2 * old.y - 1 */
  779.    else if(lold.x > ltmp.y)    /* if old.x > .5 */
  780.       lnew.x = lnew.x - ltmp.x; /* new.x = 2 * old.x - 1 */
  781.    /* end barnsley code */
  782.    return(longbailout());
  783. #endif
  784. }
  785.  
  786. SierpinskiFPFractal()
  787. {
  788.    /* following code translated from basic - see "Fractals
  789.    Everywhere" by Michael Barnsley, p. 251, Program 7.1.1 */
  790.  
  791.    new.x = old.x + old.x;
  792.    new.y = old.y + old.y;
  793.    if(old.y > .5)
  794.       new.y = new.y - 1;
  795.    else if (old.x > .5)
  796.       new.x = new.x - 1;
  797.  
  798.    /* end barnsley code */
  799.    return(floatbailout());
  800. }
  801.  
  802. LambdaexponentFractal()
  803. {
  804.    /* found this in  "Science of Fractal Images" */
  805.    FLOATEXPBAILOUT();
  806.    FPUsincos  (&old.y,&siny,&cosy);
  807.  
  808.    if (old.x >= rqlim && cosy >= 0.0) return(1);
  809.    tmpexp = exp(old.x);
  810.    tmp.x = tmpexp*cosy;
  811.    tmp.y = tmpexp*siny;
  812.  
  813.    /*multiply by lamda */
  814.    new.x = floatparm->x*tmp.x - floatparm->y*tmp.y;
  815.    new.y = floatparm->y*tmp.x + floatparm->x*tmp.y;
  816.    old = new;
  817.    return(0);
  818. }
  819.  
  820. LongLambdaexponentFractal()
  821. {
  822. #ifndef XFRACT
  823.    /* found this in  "Science of Fractal Images" */
  824.    LONGEXPBAILOUT();
  825.  
  826.    SinCos086  (lold.y, &lsiny,    &lcosy);
  827.  
  828.    if (lold.x >= llimit && lcosy >= 0L) return(1);
  829.    longtmp = Exp086(lold.x);
  830.  
  831.    ltmp.x = multiply(longtmp,       lcosy,   bitshift);
  832.    ltmp.y = multiply(longtmp,       lsiny,   bitshift);
  833.  
  834.    lnew.x  = multiply(longparm->x, ltmp.x, bitshift)
  835.        - multiply(longparm->y, ltmp.y, bitshift);
  836.    lnew.y  = multiply(longparm->x, ltmp.y, bitshift)
  837.        + multiply(longparm->y, ltmp.x, bitshift);
  838.    lold = lnew;
  839.    return(0);
  840. #endif
  841. }
  842.  
  843. FloatTrigPlusExponentFractal()
  844. {
  845.    /* another Scientific American biomorph type */
  846.    /* z(n+1) = e**z(n) + trig(z(n)) + C */
  847.  
  848.    if (fabs(old.x) >= 6.4e2) return(1); /* DOMAIN errors */
  849.    tmpexp = exp(old.x);
  850.    FPUsincos  (&old.y,&siny,&cosy);
  851.    CMPLXtrig0(old,new);
  852.  
  853.    /*new =   trig(old) + e**old + C  */
  854.    new.x += tmpexp*cosy + floatparm->x;
  855.    new.y += tmpexp*siny + floatparm->y;
  856.    return(floatbailout());
  857. }
  858.  
  859.  
  860. LongTrigPlusExponentFractal()
  861. {
  862. #ifndef XFRACT
  863.    /* calculate exp(z) */
  864.  
  865.    /* domain check for fast transcendental functions */
  866.    TRIG16CHECK(lold.x);
  867.    TRIG16CHECK(lold.y);
  868.  
  869.    longtmp = Exp086(lold.x);
  870.    SinCos086  (lold.y, &lsiny,    &lcosy);
  871.    LCMPLXtrig0(lold,lnew);
  872.    lnew.x += multiply(longtmp,      lcosy,   bitshift) + longparm->x;
  873.    lnew.y += multiply(longtmp,      lsiny,   bitshift) + longparm->y;
  874.    return(longbailout());
  875. #endif
  876. }
  877.  
  878. MarksLambdaFractal()
  879. {
  880.    /* Mark Peterson's variation of "lambda" function */
  881.  
  882.    /* Z1 = (C^(exp-1) * Z**2) + C */
  883.    ltmp.x = ltempsqrx - ltempsqry;
  884.    ltmp.y = multiply(lold.x ,lold.y ,bitshiftless1);
  885.  
  886.    lnew.x = multiply(lcoefficient.x, ltmp.x, bitshift)
  887.     - multiply(lcoefficient.y, ltmp.y, bitshift) + longparm->x;
  888.    lnew.y = multiply(lcoefficient.x, ltmp.y, bitshift)
  889.     + multiply(lcoefficient.y, ltmp.x, bitshift) + longparm->y;
  890.  
  891.    return(longbailout());
  892. }
  893.  
  894. MarksLambdafpFractal()
  895. {
  896.    /* Mark Peterson's variation of "lambda" function */
  897.  
  898.    /* Z1 = (C^(exp-1) * Z**2) + C */
  899.    tmp.x = tempsqrx - tempsqry;
  900.    tmp.y = old.x * old.y *2;
  901.  
  902.    new.x = coefficient.x * tmp.x - coefficient.y * tmp.y + floatparm->x;
  903.    new.y = coefficient.x * tmp.y + coefficient.y * tmp.x + floatparm->y;
  904.  
  905.    return(floatbailout());
  906. }
  907.  
  908.  
  909. long XXOne, FgOne, FgTwo;
  910.  
  911. UnityFractal()
  912. {
  913. #ifndef XFRACT
  914.    /* brought to you by Mark Peterson - you won't find this in any fractal
  915.       books unless they saw it here first - Mark invented it! */
  916.    XXOne = multiply(lold.x, lold.x, bitshift) + multiply(lold.y, lold.y, bitshift);
  917.    if((XXOne > FgTwo) || (labs(XXOne - FgOne) < delmin))
  918.       return(1);
  919.    lold.y = multiply(FgTwo - XXOne, lold.x, bitshift);
  920.    lold.x = multiply(FgTwo - XXOne, lold.y, bitshift);
  921.    lnew=lold;  /* TW added this line */
  922.    return(0);
  923. #endif
  924. }
  925.  
  926. #define XXOne new.x
  927.  
  928. UnityfpFractal()
  929. {
  930.    /* brought to you by Mark Peterson - you won't find this in any fractal
  931.       books unless they saw it here first - Mark invented it! */
  932.  
  933.    XXOne = sqr(old.x) + sqr(old.y);
  934.    if((XXOne > 2.0) || (fabs(XXOne - 1.0) < ddelmin))
  935.       return(1);
  936.    old.y = (2.0 - XXOne)* old.x;
  937.    old.x = (2.0 - XXOne)* old.y;
  938.    new=old;  /* TW added this line */
  939.    return(0);
  940. }
  941.  
  942. #undef XXOne
  943.  
  944. Mandel4Fractal()
  945. {
  946.    /* By writing this code, Bert has left behind the excuse "don't
  947.       know what a fractal is, just know how to make'em go fast".
  948.       Bert is hereby declared a bonafide fractal expert! Supposedly
  949.       this routine calculates the Mandelbrot/Julia set based on the
  950.       polynomial z**4 + lambda, but I wouldn't know -- can't follow
  951.       all that integer math speedup stuff - Tim */
  952.  
  953.    /* first, compute (x + iy)**2 */
  954.    lnew.x  = ltempsqrx - ltempsqry;
  955.    lnew.y = multiply(lold.x, lold.y, bitshiftless1);
  956.    if (longbailout()) return(1);
  957.  
  958.    /* then, compute ((x + iy)**2)**2 + lambda */
  959.    lnew.x  = ltempsqrx - ltempsqry + longparm->x;
  960.    lnew.y = multiply(lold.x, lold.y, bitshiftless1) + longparm->y;
  961.    return(longbailout());
  962. }
  963.  
  964. Mandel4fpFractal()
  965. {
  966.    /* first, compute (x + iy)**2 */
  967.    new.x  = tempsqrx - tempsqry;
  968.    new.y = old.x*old.y*2;
  969.    if (floatbailout()) return(1);
  970.  
  971.    /* then, compute ((x + iy)**2)**2 + lambda */
  972.    new.x  = tempsqrx - tempsqry + floatparm->x;
  973.    new.y =  old.x*old.y*2 + floatparm->y;
  974.    return(floatbailout());
  975. }
  976.  
  977. floatZtozPluszpwrFractal()
  978. {
  979.    cpower(&old,(int)param[2],&new);
  980.    old = ComplexPower(old,old);
  981.    new.x = new.x + old.x +floatparm->x;
  982.    new.y = new.y + old.y +floatparm->y;
  983.    return(floatbailout());
  984. }
  985.  
  986. longZpowerFractal()
  987. {
  988. #ifndef XFRACT
  989.    if(lcpower(&lold,c_exp,&lnew,bitshift))
  990.       lnew.x = lnew.y = 8L<<bitshift;
  991.    lnew.x += longparm->x;
  992.    lnew.y += longparm->y;
  993.    return(longbailout());
  994. #endif
  995. }
  996.  
  997. longCmplxZpowerFractal()
  998. {
  999. #ifndef XFRACT
  1000.    _CMPLX x, y;
  1001.  
  1002.    x.x = (double)lold.x / fudge;
  1003.    x.y = (double)lold.y / fudge;
  1004.    y.x = (double)lparm2.x / fudge;
  1005.    y.y = (double)lparm2.y / fudge;
  1006.    x = ComplexPower(x, y);
  1007.    if(fabs(x.x) < fgLimit && fabs(x.y) < fgLimit) {
  1008.       lnew.x = (long)(x.x * fudge);
  1009.       lnew.y = (long)(x.y * fudge);
  1010.    }
  1011.    else
  1012.       overflow = 1;
  1013.    lnew.x += longparm->x;
  1014.    lnew.y += longparm->y;
  1015.    return(longbailout());
  1016. #endif
  1017. }
  1018.  
  1019. floatZpowerFractal()
  1020. {
  1021.    cpower(&old,c_exp,&new);
  1022.    new.x += floatparm->x;
  1023.    new.y += floatparm->y;
  1024.    return(floatbailout());
  1025. }
  1026.  
  1027. floatCmplxZpowerFractal()
  1028. {
  1029.    new = ComplexPower(old, parm2);
  1030.    new.x += floatparm->x;
  1031.    new.y += floatparm->y;
  1032.    return(floatbailout());
  1033. }
  1034.  
  1035. Barnsley3Fractal()
  1036. {
  1037.    /* An unnamed Mandelbrot/Julia function from "Fractals
  1038.    Everywhere" by Michael Barnsley, p. 292, example 4.1 */
  1039.  
  1040.    /* calculate intermediate products */
  1041.    oldxinitx   = multiply(lold.x, lold.x, bitshift);
  1042.    oldyinity   = multiply(lold.y, lold.y, bitshift);
  1043.    oldxinity   = multiply(lold.x, lold.y, bitshift);
  1044.  
  1045.    /* orbit calculation */
  1046.    if(lold.x > 0)
  1047.    {
  1048.       lnew.x = oldxinitx   - oldyinity - fudge;
  1049.       lnew.y = oldxinity << 1;
  1050.    }
  1051.    else
  1052.    {
  1053.       lnew.x = oldxinitx - oldyinity - fudge
  1054.        + multiply(longparm->x,lold.x,bitshift);
  1055.       lnew.y = oldxinity <<1;
  1056.  
  1057.       /* This term added by Tim Wegner to make dependent on the
  1058.      imaginary part of the parameter. (Otherwise Mandelbrot
  1059.      is uninteresting. */
  1060.       lnew.y += multiply(longparm->y,lold.x,bitshift);
  1061.    }
  1062.    return(longbailout());
  1063. }
  1064.  
  1065. Barnsley3FPFractal()
  1066. {
  1067.    /* An unnamed Mandelbrot/Julia function from "Fractals
  1068.    Everywhere" by Michael Barnsley, p. 292, example 4.1 */
  1069.  
  1070.  
  1071.    /* calculate intermediate products */
  1072.    foldxinitx  = old.x * old.x;
  1073.    foldyinity  = old.y * old.y;
  1074.    foldxinity  = old.x * old.y;
  1075.  
  1076.    /* orbit calculation */
  1077.    if(old.x > 0)
  1078.    {
  1079.       new.x = foldxinitx - foldyinity - 1.0;
  1080.       new.y = foldxinity * 2;
  1081.    }
  1082.    else
  1083.    {
  1084.       new.x = foldxinitx - foldyinity -1.0 + floatparm->x * old.x;
  1085.       new.y = foldxinity * 2;
  1086.  
  1087.       /* This term added by Tim Wegner to make dependent on the
  1088.      imaginary part of the parameter. (Otherwise Mandelbrot
  1089.      is uninteresting. */
  1090.       new.y += floatparm->y * old.x;
  1091.    }
  1092.    return(floatbailout());
  1093. }
  1094.  
  1095. TrigPlusZsquaredFractal()
  1096. {
  1097. #ifndef XFRACT
  1098.    /* From Scientific American, July 1989 */
  1099.    /* A Biomorph              */
  1100.    /* z(n+1) = trig(z(n))+z(n)**2+C      */
  1101.    LCMPLXtrig0(lold,lnew);
  1102.    lnew.x += ltempsqrx - ltempsqry + longparm->x;
  1103.    lnew.y += multiply(lold.x, lold.y, bitshiftless1) + longparm->y;
  1104.    return(longbailout());
  1105. #endif
  1106. }
  1107.  
  1108. TrigPlusZsquaredfpFractal()
  1109. {
  1110.    /* From Scientific American, July 1989 */
  1111.    /* A Biomorph              */
  1112.    /* z(n+1) = trig(z(n))+z(n)**2+C      */
  1113.  
  1114.    CMPLXtrig0(old,new);
  1115.    new.x += tempsqrx - tempsqry + floatparm->x;
  1116.    new.y += 2.0 * old.x * old.y + floatparm->y;
  1117.    return(floatbailout());
  1118. }
  1119.  
  1120. Richard8fpFractal()
  1121. {
  1122.    /*  Richard8 {c = z = pixel: z=sin(z)+sin(pixel),|z|<=50} */
  1123.    CMPLXtrig0(old,new);
  1124. /*   CMPLXtrig1(*floatparm,tmp); */
  1125.    new.x += tmp.x;
  1126.    new.y += tmp.y;
  1127.    return(floatbailout());
  1128. }
  1129.  
  1130. Richard8Fractal()
  1131. {
  1132. #ifndef XFRACT
  1133.    /*  Richard8 {c = z = pixel: z=sin(z)+sin(pixel),|z|<=50} */
  1134.    LCMPLXtrig0(lold,lnew);
  1135. /*   LCMPLXtrig1(*longparm,ltmp); */
  1136.    lnew.x += ltmp.x;
  1137.    lnew.y += ltmp.y;
  1138.    return(longbailout());
  1139. #endif
  1140. }
  1141.  
  1142. PopcornFractal()
  1143. {
  1144.    extern int row;
  1145.    tmp = old;
  1146.    tmp.x *= 3.0;
  1147.    tmp.y *= 3.0;
  1148.    FPUsincos(&tmp.x,&sinx,&cosx);
  1149.    FPUsincos(&tmp.y,&siny,&cosy);
  1150.    tmp.x = sinx/cosx + old.x;
  1151.    tmp.y = siny/cosy + old.y;
  1152.    FPUsincos(&tmp.x,&sinx,&cosx);
  1153.    FPUsincos(&tmp.y,&siny,&cosy);
  1154.    new.x = old.x - parm.x*siny;
  1155.    new.y = old.y - parm.x*sinx;
  1156.    /*
  1157.    new.x = old.x - parm.x*sin(old.y+tan(3*old.y));
  1158.    new.y = old.y - parm.x*sin(old.x+tan(3*old.x));
  1159.    */
  1160.    if(plot == noplot)
  1161.    {
  1162.       plot_orbit(new.x,new.y,1+row%colors);
  1163.       old = new;
  1164.    }
  1165.    else
  1166.    /* FLOATBAILOUT(); */
  1167.    /* PB The above line was weird, not what it seems to be!  But, bracketing
  1168.      it or always doing it (either of which seem more likely to be what
  1169.      was intended) changes the image for the worse, so I'm not touching it.
  1170.      Same applies to int form in next routine. */
  1171.    /* PB later: recoded inline, still leaving it weird */
  1172.       tempsqrx = sqr(new.x);
  1173.    tempsqry = sqr(new.y);
  1174.    if((magnitude = tempsqrx + tempsqry) >= rqlim) return(1);
  1175.    old = new;
  1176.    return(0);
  1177. }
  1178.  
  1179. LPopcornFractal()
  1180. {
  1181. #ifndef XFRACT
  1182.    extern int row;
  1183.    ltmp = lold;
  1184.    ltmp.x *= 3L;
  1185.    ltmp.y *= 3L;
  1186.    LTRIGARG(ltmp.x);
  1187.    LTRIGARG(ltmp.y);
  1188.    SinCos086(ltmp.x,&lsinx,&lcosx);
  1189.    SinCos086(ltmp.y,&lsiny,&lcosy);
  1190.    ltmp.x = divide(lsinx,lcosx,bitshift) + lold.x;
  1191.    ltmp.y = divide(lsiny,lcosy,bitshift) + lold.y;
  1192.    LTRIGARG(ltmp.x);
  1193.    LTRIGARG(ltmp.y);
  1194.    SinCos086(ltmp.x,&lsinx,&lcosx);
  1195.    SinCos086(ltmp.y,&lsiny,&lcosy);
  1196.    lnew.x = lold.x - multiply(lparm.x,lsiny,bitshift);
  1197.    lnew.y = lold.y - multiply(lparm.x,lsinx,bitshift);
  1198.    if(plot == noplot)
  1199.    {
  1200.       iplot_orbit(lnew.x,lnew.y,1+row%colors);
  1201.       lold = lnew;
  1202.    }
  1203.    else
  1204.       LONGBAILOUT();
  1205.    /* PB above still the old way, is weird, see notes in FP popcorn case */
  1206.    return(0);
  1207. #endif
  1208. }
  1209.  
  1210. int MarksCplxMand(void)
  1211. {
  1212.    tmp.x = tempsqrx - tempsqry;
  1213.    tmp.y = 2*old.x*old.y;
  1214.    FPUcplxmul(&tmp, &Coefficient, &new);
  1215.    new.x += floatparm->x;
  1216.    new.y += floatparm->y;
  1217.    return(floatbailout());
  1218. }
  1219.  
  1220. int SpiderfpFractal(void)
  1221. {
  1222.    /* Spider(XAXIS) { c=z=pixel: z=z*z+c; c=c/2+z, |z|<=4 } */
  1223.    new.x = tempsqrx - tempsqry + tmp.x;
  1224.    new.y = 2 * old.x * old.y + tmp.y;
  1225.    tmp.x = tmp.x/2 + new.x;
  1226.    tmp.y = tmp.y/2 + new.y;
  1227.    return(floatbailout());
  1228. }
  1229.  
  1230. SpiderFractal(void)
  1231. {
  1232. #ifndef XFRACT
  1233.    /* Spider(XAXIS) { c=z=pixel: z=z*z+c; c=c/2+z, |z|<=4 } */
  1234.    lnew.x  = ltempsqrx - ltempsqry + ltmp.x;
  1235.    lnew.y = multiply(lold.x, lold.y, bitshiftless1) + ltmp.y;
  1236.    ltmp.x = (ltmp.x >> 1) + lnew.x;
  1237.    ltmp.y = (ltmp.y >> 1) + lnew.y;
  1238.    return(longbailout());
  1239. #endif
  1240. }
  1241.  
  1242. TetratefpFractal()
  1243. {
  1244.    /* Tetrate(XAXIS) { c=z=pixel: z=c^z, |z|<=(P1+3) } */
  1245.    new = ComplexPower(*floatparm,old);
  1246.    return(floatbailout());
  1247. }
  1248.  
  1249. ZXTrigPlusZFractal()
  1250. {
  1251. #ifndef XFRACT
  1252.    /* z = (p1*z*trig(z))+p2*z */
  1253.    LCMPLXtrig0(lold,ltmp);        /* ltmp  = trig(old)         */
  1254.    LCMPLXmult(lparm,ltmp,ltmp);      /* ltmp  = p1*trig(old)          */
  1255.    LCMPLXmult(lold,ltmp,ltmp2);      /* ltmp2 = p1*old*trig(old)      */
  1256.    LCMPLXmult(lparm2,lold,ltmp);     /* ltmp  = p2*old              */
  1257.    LCMPLXadd(ltmp2,ltmp,lnew);         /* lnew  = p1*trig(old) + p2*old */
  1258.    return(longbailout());
  1259. #endif
  1260. }
  1261.  
  1262. ScottZXTrigPlusZFractal()
  1263. {
  1264. #ifndef XFRACT
  1265.    /* z = (z*trig(z))+z */
  1266.    LCMPLXtrig0(lold,ltmp);        /* ltmp  = trig(old)       */
  1267.    LCMPLXmult(lold,ltmp,lnew);         /* lnew  = old*trig(old)    */
  1268.    LCMPLXadd(lnew,lold,lnew);         /* lnew  = trig(old) + old */
  1269.    return(longbailout());
  1270. #endif
  1271. }
  1272.  
  1273. SkinnerZXTrigSubZFractal()
  1274. {
  1275. #ifndef XFRACT
  1276.    /* z = (z*trig(z))-z */
  1277.    LCMPLXtrig0(lold,ltmp);        /* ltmp  = trig(old)       */
  1278.    LCMPLXmult(lold,ltmp,lnew);         /* lnew  = old*trig(old)    */
  1279.    LCMPLXsub(lnew,lold,lnew);         /* lnew  = trig(old) - old */
  1280.    return(longbailout());
  1281. #endif
  1282. }
  1283.  
  1284. ZXTrigPlusZfpFractal()
  1285. {
  1286.    /* z = (p1*z*trig(z))+p2*z */
  1287.    CMPLXtrig0(old,tmp);      /* tmp  = trig(old)         */
  1288.    CMPLXmult(parm,tmp,tmp);     /* tmp  = p1*trig(old)      */
  1289.    CMPLXmult(old,tmp,tmp2);     /* tmp2 = p1*old*trig(old)     */
  1290.    CMPLXmult(parm2,old,tmp);     /* tmp  = p2*old         */
  1291.    CMPLXadd(tmp2,tmp,new);     /* new  = p1*trig(old) + p2*old */
  1292.    return(floatbailout());
  1293. }
  1294.  
  1295. ScottZXTrigPlusZfpFractal()
  1296. {
  1297.    /* z = (z*trig(z))+z */
  1298.    CMPLXtrig0(old,tmp);     /* tmp    = trig(old)      */
  1299.    CMPLXmult(old,tmp,new);     /* new  = old*trig(old)   */
  1300.    CMPLXadd(new,old,new);     /* new  = trig(old) + old */
  1301.    return(floatbailout());
  1302. }
  1303.  
  1304. SkinnerZXTrigSubZfpFractal()
  1305. {
  1306.    /* z = (z*trig(z))-z */
  1307.    CMPLXtrig0(old,tmp);     /* tmp    = trig(old)      */
  1308.    CMPLXmult(old,tmp,new);     /* new  = old*trig(old)   */
  1309.    CMPLXsub(new,old,new);     /* new  = trig(old) - old */
  1310.    return(floatbailout());
  1311. }
  1312.  
  1313. Sqr1overTrigFractal()
  1314. {
  1315. #ifndef XFRACT
  1316.    /* z = sqr(1/trig(z)) */
  1317.    LCMPLXtrig0(lold,lold);
  1318.    LCMPLXrecip(lold,lold);
  1319.    LCMPLXsqr(lold,lnew);
  1320.    return(longbailout());
  1321. #endif
  1322. }
  1323.  
  1324. Sqr1overTrigfpFractal()
  1325. {
  1326.    /* z = sqr(1/trig(z)) */
  1327.    CMPLXtrig0(old,old);
  1328.    CMPLXrecip(old,old);
  1329.    CMPLXsqr(old,new);
  1330.    return(floatbailout());
  1331. }
  1332.  
  1333. TrigPlusTrigFractal()
  1334. {
  1335. #ifndef XFRACT
  1336.    /* z = trig(0,z)*p1+trig1(z)*p2 */
  1337.    LCMPLXtrig0(lold,ltmp);
  1338.    LCMPLXmult(lparm,ltmp,ltmp);
  1339.    LCMPLXtrig1(lold,ltmp2);
  1340.    LCMPLXmult(lparm2,ltmp2,lold);
  1341.    LCMPLXadd(ltmp,lold,lnew);
  1342.    return(longbailout());
  1343. #endif
  1344. }
  1345.  
  1346. TrigPlusTrigfpFractal()
  1347. {
  1348.    /* z = trig0(z)*p1+trig1(z)*p2 */
  1349.    CMPLXtrig0(old,tmp);
  1350.    CMPLXmult(parm,tmp,tmp);
  1351.    CMPLXtrig1(old,old);
  1352.    CMPLXmult(parm2,old,old);
  1353.    CMPLXadd(tmp,old,new);
  1354.    return(floatbailout());
  1355. }
  1356.  
  1357. /* The following four fractals are based on the idea of parallel
  1358.    or alternate calculations.  The shift is made when the mod
  1359.    reaches a given value.  JCO  5/6/92 */
  1360.  
  1361. LambdaTrigOrTrigFractal()
  1362. {
  1363. #ifndef XFRACT
  1364.    /* z = trig0(z)*p1 if mod(old) < p2.x and
  1365.           trig1(z)*p1 if mod(old) >= p2.x */
  1366.    if ((LCMPLXmod(lold)) < lparm2.x){
  1367.      LCMPLXtrig0(lold,ltmp);
  1368.      LCMPLXmult(*longparm,ltmp,lnew);}
  1369.    else{
  1370.      LCMPLXtrig1(lold,ltmp);
  1371.      LCMPLXmult(*longparm,ltmp,lnew);}
  1372.    return(longbailout());
  1373. #endif
  1374. }
  1375.  
  1376. LambdaTrigOrTrigfpFractal()
  1377. {
  1378.    /* z = trig0(z)*p1 if mod(old) < p2.x and
  1379.           trig1(z)*p1 if mod(old) >= p2.x */
  1380.    if (CMPLXmod(old) < parm2.x){
  1381.      CMPLXtrig0(old,old);
  1382.      FPUcplxmul(floatparm,&old,&new);}
  1383.    else{
  1384.      CMPLXtrig1(old,old);
  1385.      FPUcplxmul(floatparm,&old,&new);}
  1386.    return(floatbailout());
  1387. }
  1388.  
  1389. JuliaTrigOrTrigFractal()
  1390. {
  1391. #ifndef XFRACT
  1392.    /* z = trig0(z)+p1 if mod(old) < p2.x and
  1393.           trig1(z)+p1 if mod(old) >= p2.x */
  1394.    if (LCMPLXmod(lold) < lparm2.x){
  1395.      LCMPLXtrig0(lold,ltmp);
  1396.      LCMPLXadd(*longparm,ltmp,lnew);}
  1397.    else{
  1398.      LCMPLXtrig1(lold,ltmp);
  1399.      LCMPLXadd(*longparm,ltmp,lnew);}
  1400.    return(longbailout());
  1401. #endif
  1402. }
  1403.  
  1404. JuliaTrigOrTrigfpFractal()
  1405. {
  1406.    /* z = trig0(z)+p1 if mod(old) < p2.x and
  1407.           trig1(z)+p1 if mod(old) >= p2.x */
  1408.    if (CMPLXmod(old) < parm2.x){
  1409.      CMPLXtrig0(old,old);
  1410.      CMPLXadd(*floatparm,old,new);}
  1411.    else{
  1412.      CMPLXtrig1(old,old);
  1413.      CMPLXadd(*floatparm,old,new);}
  1414.    return(floatbailout());
  1415. }
  1416.  
  1417. static int AplusOne, Ap1deg;
  1418. static struct MP mpAplusOne, mpAp1deg, mptmpparmy;
  1419.  
  1420. int MPCHalleyFractal()
  1421. {
  1422. #ifndef XFRACT
  1423.    /*  X(X^a - 1) = 0, Halley Map */
  1424.    /*  a = parm.x,  relaxation coeff. = parm.y,  epsilon = parm2.x  */
  1425.  
  1426. int ihal;
  1427. struct MPC mpcXtoAlessOne, mpcXtoA;
  1428. struct MPC mpcXtoAplusOne; /* a-1, a, a+1 */
  1429. struct MPC mpcFX, mpcF1prime, mpcF2prime, mpcHalnumer1;
  1430. struct MPC mpcHalnumer2, mpcHaldenom, mpctmp;
  1431.  
  1432.    MPOverflow = 0;
  1433.    mpcXtoAlessOne.x = mpcold.x;
  1434.    mpcXtoAlessOne.y = mpcold.y;
  1435.    for(ihal=2; ihal<degree; ihal++) {
  1436.      mpctmp.x = *pMPsub(*pMPmul(mpcXtoAlessOne.x,mpcold.x),*pMPmul(mpcXtoAlessOne.y,mpcold.y));
  1437.      mpctmp.y = *pMPadd(*pMPmul(mpcXtoAlessOne.x,mpcold.y),*pMPmul(mpcXtoAlessOne.y,mpcold.x));
  1438.      mpcXtoAlessOne.x = mpctmp.x;
  1439.      mpcXtoAlessOne.y = mpctmp.y;
  1440.    }
  1441.    mpcXtoA.x = *pMPsub(*pMPmul(mpcXtoAlessOne.x,mpcold.x),*pMPmul(mpcXtoAlessOne.y,mpcold.y));
  1442.    mpcXtoA.y = *pMPadd(*pMPmul(mpcXtoAlessOne.x,mpcold.y),*pMPmul(mpcXtoAlessOne.y,mpcold.x));
  1443.    mpcXtoAplusOne.x = *pMPsub(*pMPmul(mpcXtoA.x,mpcold.x),*pMPmul(mpcXtoA.y,mpcold.y));
  1444.    mpcXtoAplusOne.y = *pMPadd(*pMPmul(mpcXtoA.x,mpcold.y),*pMPmul(mpcXtoA.y,mpcold.x));
  1445.  
  1446.    mpcFX.x = *pMPsub(mpcXtoAplusOne.x, mpcold.x);
  1447.    mpcFX.y = *pMPsub(mpcXtoAplusOne.y, mpcold.y); /* FX = X^(a+1) - X  = F */
  1448.  
  1449.    mpcF2prime.x = *pMPmul(mpAp1deg, mpcXtoAlessOne.x); /* mpAp1deg in setup */
  1450.    mpcF2prime.y = *pMPmul(mpAp1deg, mpcXtoAlessOne.y);        /* F" */
  1451.  
  1452.    mpcF1prime.x = *pMPsub(*pMPmul(mpAplusOne, mpcXtoA.x), mpone);
  1453.    mpcF1prime.y = *pMPmul(mpAplusOne, mpcXtoA.y);                   /*  F'  */
  1454.  
  1455.    mpctmp.x = *pMPsub(*pMPmul(mpcF2prime.x,mpcFX.x),*pMPmul(mpcF2prime.y,mpcFX.y));
  1456.    mpctmp.y = *pMPadd(*pMPmul(mpcF2prime.x,mpcFX.y),*pMPmul(mpcF2prime.y,mpcFX.x));
  1457.    /*  F * F"  */
  1458.  
  1459.    mpcHaldenom.x = *pMPadd(mpcF1prime.x, mpcF1prime.x);
  1460.    mpcHaldenom.y = *pMPadd(mpcF1prime.y, mpcF1prime.y);      /*  2 * F'  */
  1461.  
  1462.    mpcHalnumer1 = MPCdiv(mpctmp, mpcHaldenom);        /*  F"F/2F'  */
  1463.    mpctmp.x = *pMPsub(mpcF1prime.x, mpcHalnumer1.x);
  1464.    mpctmp.y = *pMPsub(mpcF1prime.y, mpcHalnumer1.y); /*  F' - F"F/2F'  */
  1465.    mpcHalnumer2 = MPCdiv(mpcFX, mpctmp);
  1466.  
  1467.    mpctmp.x = *pMPmul(mptmpparmy,mpcHalnumer2.x); /* mptmpparmy is */
  1468.    mpctmp.y = *pMPmul(mptmpparmy,mpcHalnumer2.y); /* relaxation coef. */
  1469.    mpcnew.x = *pMPsub(mpcold.x, mpctmp.x);
  1470.    mpcnew.y = *pMPsub(mpcold.y, mpctmp.y);
  1471.  
  1472.    new.x = *pMP2d(mpcnew.x);
  1473.    new.y = *pMP2d(mpcnew.y);
  1474.  
  1475.    return(MPCHalleybailout()||MPOverflow);
  1476. #endif
  1477. }
  1478.  
  1479. HalleyFractal()
  1480. {
  1481.    /*  X(X^a - 1) = 0, Halley Map */
  1482.    /*  a = parm.x = degree, relaxation coeff. = parm.y, epsilon = parm2.x  */
  1483.  
  1484. int ihal;
  1485. _CMPLX XtoAlessOne, XtoA, XtoAplusOne; /* a-1, a, a+1 */
  1486. _CMPLX FX, F1prime, F2prime, Halnumer1, Halnumer2, Haldenom;
  1487.  
  1488.    XtoAlessOne = old;
  1489.    for(ihal=2; ihal<degree; ihal++) {
  1490.      FPUcplxmul(&old, &XtoAlessOne, &XtoAlessOne);
  1491.    }
  1492.    FPUcplxmul(&old, &XtoAlessOne, &XtoA);
  1493.    FPUcplxmul(&old, &XtoA, &XtoAplusOne);
  1494.  
  1495.    CMPLXsub(XtoAplusOne, old, FX);        /* FX = X^(a+1) - X  = F */
  1496.    F2prime.x = Ap1deg * XtoAlessOne.x; /* Ap1deg in setup */
  1497.    F2prime.y = Ap1deg * XtoAlessOne.y;        /* F" */
  1498.  
  1499.    F1prime.x = AplusOne * XtoA.x - 1.0;
  1500.    F1prime.y = AplusOne * XtoA.y;                             /*  F'  */
  1501.  
  1502.    FPUcplxmul(&F2prime, &FX, &Halnumer1);                  /*  F * F"  */
  1503.    Haldenom.x = F1prime.x + F1prime.x;
  1504.    Haldenom.y = F1prime.y + F1prime.y;                     /*  2 * F'  */
  1505.  
  1506.    FPUcplxdiv(&Halnumer1, &Haldenom, &Halnumer1);         /*  F"F/2F'  */
  1507.    CMPLXsub(F1prime, Halnumer1, Halnumer2);          /*  F' - F"F/2F'  */
  1508.    FPUcplxdiv(&FX, &Halnumer2, &Halnumer2);
  1509.    new.x = old.x - (parm.y * Halnumer2.x); /* parm.y is relaxation coef. */
  1510.    new.y = old.y - (parm.y * Halnumer2.y);
  1511.  
  1512.    return(Halleybailout());
  1513. }
  1514.  
  1515. LongPhoenixFractal()
  1516. {
  1517. #ifndef XFRACT
  1518. /* z(n+1) = z(n)^2 + p + qy(n),  y(n+1) = z(n) */
  1519.    ltmp.x = multiply(lold.x, lold.y, bitshift);
  1520.    lnew.x = ltempsqrx-ltempsqry+longparm->x+multiply(longparm->y,ltmp2.x,bitshift);
  1521.    lnew.y = (ltmp.x + ltmp.x) + multiply(longparm->y,ltmp2.y,bitshift);
  1522.    ltmp2 = lold; /* set ltmp2 to Y value */
  1523.    return(longbailout());
  1524. #endif
  1525. }
  1526.  
  1527. PhoenixFractal()
  1528. {
  1529. /* z(n+1) = z(n)^2 + p + qy(n),  y(n+1) = z(n) */
  1530.    tmp.x = old.x * old.y;
  1531.    new.x = tempsqrx - tempsqry + floatparm->x + (floatparm->y * tmp2.x);
  1532.    new.y = (tmp.x + tmp.x) + (floatparm->y * tmp2.y);
  1533.    tmp2 = old; /* set tmp2 to Y value */
  1534.    return(floatbailout());
  1535. }
  1536.  
  1537. LongPhoenixPlusFractal()
  1538. {
  1539. #ifndef XFRACT
  1540. /* z(n+1) = z(n)^(degree-1) * (z(n) + p) + qy(n),  y(n+1) = z(n) */
  1541. int i;
  1542. _LCMPLX loldplus, lnewminus;
  1543.    loldplus = lold;
  1544.    ltmp = lold;
  1545.    for(i=1; i<degree; i++) { /* degree >= 2, degree=degree-1 in setup */
  1546.       LCMPLXmult(lold,ltmp,ltmp); /* = old^(degree-1) */
  1547.    }
  1548.    loldplus.x += longparm->x;
  1549.    LCMPLXmult(ltmp, loldplus, lnewminus);
  1550.    lnew.x = lnewminus.x + multiply(longparm->y,ltmp2.x,bitshift);
  1551.    lnew.y = lnewminus.y + multiply(longparm->y,ltmp2.y,bitshift);
  1552.    ltmp2 = lold; /* set ltmp2 to Y value */
  1553.    return(longbailout());
  1554. #endif
  1555. }
  1556.  
  1557. PhoenixPlusFractal()
  1558. {
  1559. /* z(n+1) = z(n)^(degree-1) * (z(n) + p) + qy(n),  y(n+1) = z(n) */
  1560. int i;
  1561. _CMPLX oldplus, newminus;
  1562.    oldplus = old;
  1563.    tmp = old;
  1564.    for(i=1; i<degree; i++) { /* degree >= 2, degree=degree-1 in setup */
  1565.      FPUcplxmul(&old, &tmp, &tmp); /* = old^(degree-1) */
  1566.    }
  1567.    oldplus.x += floatparm->x;
  1568.    FPUcplxmul(&tmp, &oldplus, &newminus);
  1569.    new.x = newminus.x + (floatparm->y * tmp2.x);
  1570.    new.y = newminus.y + (floatparm->y * tmp2.y);
  1571.    tmp2 = old; /* set tmp2 to Y value */
  1572.    return(floatbailout());
  1573. }
  1574.  
  1575. LongPhoenixMinusFractal()
  1576. {
  1577. #ifndef XFRACT
  1578. /* z(n+1) = z(n)^(degree-2) * (z(n)^2 + p) + qy(n),  y(n+1) = z(n) */
  1579. int i;
  1580. _LCMPLX loldsqr, lnewminus;
  1581.    LCMPLXmult(lold,lold,loldsqr);
  1582.    ltmp = lold;
  1583.    for(i=1; i<degree; i++) { /* degree >= 3, degree=degree-2 in setup */
  1584.       LCMPLXmult(lold,ltmp,ltmp); /* = old^(degree-2) */
  1585.    }
  1586.    loldsqr.x += longparm->x;
  1587.    LCMPLXmult(ltmp, loldsqr, lnewminus);
  1588.    lnew.x = lnewminus.x + multiply(longparm->y,ltmp2.x,bitshift);
  1589.    lnew.y = lnewminus.y + multiply(longparm->y,ltmp2.y,bitshift);
  1590.    ltmp2 = lold; /* set ltmp2 to Y value */
  1591.    return(longbailout());
  1592. #endif
  1593. }
  1594.  
  1595. PhoenixMinusFractal()
  1596. {
  1597. /* z(n+1) = z(n)^(degree-2) * (z(n)^2 + p) + qy(n),  y(n+1) = z(n) */
  1598. int i;
  1599. _CMPLX oldsqr, newminus;
  1600.    FPUcplxmul(&old, &old, &oldsqr);
  1601.    tmp = old;
  1602.    for(i=1; i<degree; i++) { /* degree >= 3, degree=degree-2 in setup */
  1603.      FPUcplxmul(&old, &tmp, &tmp); /* = old^(degree-2) */
  1604.    }
  1605.    oldsqr.x += floatparm->x;
  1606.    FPUcplxmul(&tmp, &oldsqr, &newminus);
  1607.    new.x = newminus.x + (floatparm->y * tmp2.x);
  1608.    new.y = newminus.y + (floatparm->y * tmp2.y);
  1609.    tmp2 = old; /* set tmp2 to Y value */
  1610.    return(floatbailout());
  1611. }
  1612.  
  1613. ScottTrigPlusTrigFractal()
  1614. {
  1615. #ifndef XFRACT
  1616.    /* z = trig0(z)+trig1(z) */
  1617.    LCMPLXtrig0(lold,ltmp);
  1618.    LCMPLXtrig1(lold,lold);
  1619.    LCMPLXadd(ltmp,lold,lnew);
  1620.    return(longbailout());
  1621. #endif
  1622. }
  1623.  
  1624. ScottTrigPlusTrigfpFractal()
  1625. {
  1626.    /* z = trig0(z)+trig1(z) */
  1627.    CMPLXtrig0(old,tmp);
  1628.    CMPLXtrig1(old,tmp2);
  1629.    CMPLXadd(tmp,tmp2,new);
  1630.    return(floatbailout());
  1631. }
  1632.  
  1633. SkinnerTrigSubTrigFractal()
  1634. {
  1635. #ifndef XFRACT
  1636.    /* z = trig(0,z)-trig1(z) */
  1637.    LCMPLXtrig0(lold,ltmp);
  1638.    LCMPLXtrig1(lold,ltmp2);
  1639.    LCMPLXsub(ltmp,ltmp2,lnew);
  1640.    return(longbailout());
  1641. #endif
  1642. }
  1643.  
  1644. SkinnerTrigSubTrigfpFractal()
  1645. {
  1646.    /* z = trig0(z)-trig1(z) */
  1647.    CMPLXtrig0(old,tmp);
  1648.    CMPLXtrig1(old,tmp2);
  1649.    CMPLXsub(tmp,tmp2,new);
  1650.    return(floatbailout());
  1651. }
  1652.  
  1653.  
  1654. TrigXTrigfpFractal()
  1655. {
  1656.    /* z = trig0(z)*trig1(z) */
  1657.    CMPLXtrig0(old,tmp);
  1658.    CMPLXtrig1(old,old);
  1659.    CMPLXmult(tmp,old,new);
  1660.    return(floatbailout());
  1661. }
  1662.  
  1663. TrigXTrigFractal()
  1664. {
  1665. #ifndef XFRACT
  1666.    _LCMPLX ltmp2;
  1667.    /* z = trig0(z)*trig1(z) */
  1668.    LCMPLXtrig0(lold,ltmp);
  1669.    LCMPLXtrig1(lold,ltmp2);
  1670.    LCMPLXmult(ltmp,ltmp2,lnew);
  1671.    if(overflow)
  1672.       TryFloatFractal(TrigXTrigfpFractal);
  1673.    return(longbailout());
  1674. #endif
  1675. }
  1676.  
  1677.  /* call float version of fractal if integer math overflow */
  1678. TryFloatFractal(int (*fpFractal)())
  1679. {
  1680.    overflow=0;
  1681.    /* lold had better not be changed! */
  1682.    old.x = lold.x; old.x /= fudge;
  1683.    old.y = lold.y; old.y /= fudge;
  1684.    tempsqrx = sqr(old.x);
  1685.    tempsqry = sqr(old.y);
  1686.    fpFractal();
  1687.    lnew.x = new.x/fudge;
  1688.    lnew.y = new.y/fudge;
  1689.    return(0);
  1690. }
  1691.  
  1692. /********************************************************************/
  1693. /*  Next six orbit functions are one type - extra functions are     */
  1694. /*    special cases written for speed.                    */
  1695. /********************************************************************/
  1696.  
  1697. TrigPlusSqrFractal() /* generalization of Scott and Skinner types */
  1698. {
  1699. #ifndef XFRACT
  1700.    /* { z=pixel: z=(p1,p2)*trig(z)+(p3,p4)*sqr(z), |z|<BAILOUT } */
  1701.    LCMPLXtrig0(lold,ltmp);     /* ltmp = trig(lold)               */
  1702.    LCMPLXmult(lparm,ltmp,lnew); /* lnew = lparm*trig(lold)            */
  1703.    LCMPLXsqr_old(ltmp);     /* ltmp = sqr(lold)                */
  1704.    LCMPLXmult(lparm2,ltmp,ltmp);/* ltmp = lparm2*sqr(lold)            */
  1705.    LCMPLXadd(lnew,ltmp,lnew);    /* lnew = lparm*trig(lold)+lparm2*sqr(lold) */
  1706.    return(longbailout());
  1707. #endif
  1708. }
  1709.  
  1710. TrigPlusSqrfpFractal() /* generalization of Scott and Skinner types */
  1711. {
  1712.    /* { z=pixel: z=(p1,p2)*trig(z)+(p3,p4)*sqr(z), |z|<BAILOUT } */
  1713.    CMPLXtrig0(old,tmp);     /* tmp = trig(old)               */
  1714.    CMPLXmult(parm,tmp,new); /* new = parm*trig(old)           */
  1715.    CMPLXsqr_old(tmp);         /* tmp = sqr(old)                */
  1716.    CMPLXmult(parm2,tmp,tmp2); /* tmp = parm2*sqr(old)             */
  1717.    CMPLXadd(new,tmp2,new);    /* new = parm*trig(old)+parm2*sqr(old) */
  1718.    return(floatbailout());
  1719. }
  1720.  
  1721. ScottTrigPlusSqrFractal()
  1722. {
  1723. #ifndef XFRACT
  1724.    /*  { z=pixel: z=trig(z)+sqr(z), |z|<BAILOUT } */
  1725.    LCMPLXtrig0(lold,lnew);    /* lnew = trig(lold)         */
  1726.    LCMPLXsqr_old(ltmp);        /* lold = sqr(lold)          */
  1727.    LCMPLXadd(ltmp,lnew,lnew);  /* lnew = trig(lold)+sqr(lold) */
  1728.    return(longbailout());
  1729. #endif
  1730. }
  1731.  
  1732. ScottTrigPlusSqrfpFractal() /* float version */
  1733. {
  1734.    /* { z=pixel: z=sin(z)+sqr(z), |z|<BAILOUT } */
  1735.    CMPLXtrig0(old,new);       /* new = trig(old)      */
  1736.    CMPLXsqr_old(tmp);           /* tmp = sqr(old)       */
  1737.    CMPLXadd(new,tmp,new);      /* new = trig(old)+sqr(old) */
  1738.    return(floatbailout());
  1739. }
  1740.  
  1741. SkinnerTrigSubSqrFractal()
  1742. {
  1743. #ifndef XFRACT
  1744.    /* { z=pixel: z=sin(z)-sqr(z), |z|<BAILOUT }           */
  1745.    LCMPLXtrig0(lold,lnew);    /* lnew = trig(lold)         */
  1746.    LCMPLXsqr_old(ltmp);        /* lold = sqr(lold)          */
  1747.    LCMPLXsub(lnew,ltmp,lnew);  /* lnew = trig(lold)-sqr(lold) */
  1748.    return(longbailout());
  1749. #endif
  1750. }
  1751.  
  1752. SkinnerTrigSubSqrfpFractal()
  1753. {
  1754.    /* { z=pixel: z=sin(z)-sqr(z), |z|<BAILOUT } */
  1755.    CMPLXtrig0(old,new);       /* new = trig(old) */
  1756.    CMPLXsqr_old(tmp);           /* old = sqr(old)  */
  1757.    CMPLXsub(new,tmp,new);      /* new = trig(old)-sqr(old) */
  1758.    return(floatbailout());
  1759. }
  1760.  
  1761. TrigZsqrdfpFractal()
  1762. {
  1763.    /* { z=pixel: z=trig(z*z), |z|<TEST } */
  1764.    CMPLXsqr_old(tmp);
  1765.    CMPLXtrig0(tmp,new);
  1766.    return(floatbailout());
  1767. }
  1768.  
  1769. TrigZsqrdFractal() /* this doesn't work very well */
  1770. {
  1771. #ifndef XFRACT
  1772.    /* { z=pixel: z=trig(z*z), |z|<TEST } */
  1773.    LCMPLXsqr_old(ltmp);
  1774.    LCMPLXtrig0(ltmp,lnew);
  1775.    if(overflow)
  1776.       TryFloatFractal(TrigZsqrdfpFractal);
  1777.    return(longbailout());
  1778. #endif
  1779. }
  1780.  
  1781. SqrTrigFractal()
  1782. {
  1783. #ifndef XFRACT
  1784.    /* { z=pixel: z=sqr(trig(z)), |z|<TEST} */
  1785.    LCMPLXtrig0(lold,ltmp);
  1786.    LCMPLXsqr(ltmp,lnew);
  1787.    return(longbailout());
  1788. #endif
  1789. }
  1790.  
  1791. SqrTrigfpFractal()
  1792. {
  1793.    /* SZSB(XYAXIS) { z=pixel, TEST=(p1+3): z=sin(z)*sin(z), |z|<TEST} */
  1794.    CMPLXtrig0(old,tmp);
  1795.    CMPLXsqr(tmp,new);
  1796.    return(floatbailout());
  1797. }
  1798.  
  1799.  
  1800. Magnet1Fractal()    /*      Z = ((Z**2 + C - 1)/(2Z + C - 2))**2      */
  1801.   {              /*  In "Beauty of Fractals", code by Kev Allen. */
  1802.     _CMPLX top, bot, tmp;
  1803.     double div;
  1804.  
  1805.     top.x = tempsqrx - tempsqry + floatparm->x - 1; /* top = Z**2+C-1 */
  1806.     top.y = old.x * old.y;
  1807.     top.y = top.y + top.y + floatparm->y;
  1808.  
  1809.     bot.x = old.x + old.x + floatparm->x - 2;        /* bot = 2*Z+C-2  */
  1810.     bot.y = old.y + old.y + floatparm->y;
  1811.  
  1812.     div = bot.x*bot.x + bot.y*bot.y;            /* tmp = top/bot  */
  1813.     if (div < FLT_MIN) return(1);
  1814.     tmp.x = (top.x*bot.x + top.y*bot.y)/div;
  1815.     tmp.y = (top.y*bot.x - top.x*bot.y)/div;
  1816.  
  1817.     new.x = (tmp.x + tmp.y) * (tmp.x - tmp.y);        /* Z = tmp**2     */
  1818.     new.y = tmp.x * tmp.y;
  1819.     new.y += new.y;
  1820.  
  1821.     return(floatbailout());
  1822.   }
  1823.  
  1824. Magnet2Fractal()  /* Z = ((Z**3 + 3(C-1)Z + (C-1)(C-2)    ) /     */
  1825.             /*         (3Z**2 + 3(C-2)Z + (C-1)(C-2)+1) )**2  */
  1826.   {            /*     In "Beauty of Fractals", code by Kev Allen.  */
  1827.     _CMPLX top, bot, tmp;
  1828.     double div;
  1829.  
  1830.     top.x = old.x * (tempsqrx-tempsqry-tempsqry-tempsqry + T_Cm1.x)
  1831.       - old.y * T_Cm1.y + T_Cm1Cm2.x;
  1832.     top.y = old.y * (tempsqrx+tempsqrx+tempsqrx-tempsqry + T_Cm1.x)
  1833.       + old.x * T_Cm1.y + T_Cm1Cm2.y;
  1834.  
  1835.     bot.x = tempsqrx - tempsqry;
  1836.     bot.x = bot.x + bot.x + bot.x
  1837.       + old.x * T_Cm2.x - old.y * T_Cm2.y
  1838.       + T_Cm1Cm2.x + 1.0;
  1839.     bot.y = old.x * old.y;
  1840.     bot.y += bot.y;
  1841.     bot.y = bot.y + bot.y + bot.y
  1842.       + old.x * T_Cm2.y + old.y * T_Cm2.x
  1843.       + T_Cm1Cm2.y;
  1844.  
  1845.     div = bot.x*bot.x + bot.y*bot.y;            /* tmp = top/bot  */
  1846.     if (div < FLT_MIN) return(1);
  1847.     tmp.x = (top.x*bot.x + top.y*bot.y)/div;
  1848.     tmp.y = (top.y*bot.x - top.x*bot.y)/div;
  1849.  
  1850.     new.x = (tmp.x + tmp.y) * (tmp.x - tmp.y);        /* Z = tmp**2     */
  1851.     new.y = tmp.x * tmp.y;
  1852.     new.y += new.y;
  1853.  
  1854.     return(floatbailout());
  1855.   }
  1856.  
  1857. LambdaTrigFractal()
  1858. {
  1859. #ifndef XFRACT
  1860.    LONGXYTRIGBAILOUT();
  1861.    LCMPLXtrig0(lold,ltmp);         /* ltmp = trig(lold)        */
  1862.    LCMPLXmult(*longparm,ltmp,lnew);   /* lnew = longparm*trig(lold)  */
  1863.    lold = lnew;
  1864.    return(0);
  1865. #endif
  1866. }
  1867.  
  1868. LambdaTrigfpFractal()
  1869. {
  1870.    FLOATXYTRIGBAILOUT();
  1871.    CMPLXtrig0(old,tmp);          /* tmp = trig(old)       */
  1872.    CMPLXmult(*floatparm,tmp,new);   /* new = longparm*trig(old)  */
  1873.    old = new;
  1874.    return(0);
  1875. }
  1876.  
  1877. /* bailouts are different for different trig functions */
  1878. LambdaTrigFractal1()
  1879. {
  1880. #ifndef XFRACT
  1881.    LONGTRIGBAILOUT(); /* sin,cos */
  1882.    LCMPLXtrig0(lold,ltmp);         /* ltmp = trig(lold)        */
  1883.    LCMPLXmult(*longparm,ltmp,lnew);   /* lnew = longparm*trig(lold)  */
  1884.    lold = lnew;
  1885.    return(0);
  1886. #endif
  1887. }
  1888.  
  1889. LambdaTrigfpFractal1()
  1890. {
  1891.    FLOATTRIGBAILOUT(); /* sin,cos */
  1892.    CMPLXtrig0(old,tmp);          /* tmp = trig(old)       */
  1893.    CMPLXmult(*floatparm,tmp,new);   /* new = longparm*trig(old)  */
  1894.    old = new;
  1895.    return(0);
  1896. }
  1897.  
  1898. LambdaTrigFractal2()
  1899. {
  1900. #ifndef XFRACT
  1901.    LONGHTRIGBAILOUT(); /* sinh,cosh */
  1902.    LCMPLXtrig0(lold,ltmp);         /* ltmp = trig(lold)        */
  1903.    LCMPLXmult(*longparm,ltmp,lnew);   /* lnew = longparm*trig(lold)  */
  1904.    lold = lnew;
  1905.    return(0);
  1906. #endif
  1907. }
  1908.  
  1909. LambdaTrigfpFractal2()
  1910. {
  1911. #ifndef XFRACT
  1912.    FLOATHTRIGBAILOUT(); /* sinh,cosh */
  1913.    CMPLXtrig0(old,tmp);          /* tmp = trig(old)       */
  1914.    CMPLXmult(*floatparm,tmp,new);   /* new = longparm*trig(old)  */
  1915.    old = new;
  1916.    return(0);
  1917. #endif
  1918. }
  1919.  
  1920. ManOWarFractal()
  1921. {
  1922. #ifndef XFRACT
  1923.    /* From Art Matrix via Lee Skinner */
  1924.    lnew.x  = ltempsqrx - ltempsqry + ltmp.x + longparm->x;
  1925.    lnew.y = multiply(lold.x, lold.y, bitshiftless1) + ltmp.y + longparm->y;
  1926.    ltmp = lold;
  1927.    return(longbailout());
  1928. #endif
  1929. }
  1930.  
  1931. ManOWarfpFractal()
  1932. {
  1933.    /* From Art Matrix via Lee Skinner */
  1934.    /* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
  1935.    new.x = tempsqrx - tempsqry + tmp.x + floatparm->x;
  1936.    new.y = 2.0 * old.x * old.y + tmp.y + floatparm->y;
  1937.    tmp = old;
  1938.    return(floatbailout());
  1939. }
  1940.  
  1941. /*
  1942.    MarksMandelPwr (XAXIS) {
  1943.       z = pixel, c = z ^ (z - 1):
  1944.      z = c * sqr(z) + pixel,
  1945.       |z| <= 4
  1946.    }
  1947. */
  1948.  
  1949. MarksMandelPwrfpFractal()
  1950. {
  1951.    CMPLXtrig0(old,new);
  1952.    CMPLXmult(tmp,new,new);
  1953.    new.x += floatparm->x;
  1954.    new.y += floatparm->y;
  1955.    return(floatbailout());
  1956. }
  1957.  
  1958. MarksMandelPwrFractal()
  1959. {
  1960. #ifndef XFRACT
  1961.    LCMPLXtrig0(lold,lnew);
  1962.    LCMPLXmult(ltmp,lnew,lnew);
  1963.    lnew.x += longparm->x;
  1964.    lnew.y += longparm->y;
  1965.    return(longbailout());
  1966. #endif
  1967. }
  1968.  
  1969. /* I was coding Marksmandelpower and failed to use some temporary
  1970.    variables. The result was nice, and since my name is not on any fractal,
  1971.    I thought I would immortalize myself with this error!
  1972.         Tim Wegner */
  1973.  
  1974.  
  1975. TimsErrorfpFractal()
  1976. {
  1977.    CMPLXtrig0(old,new);
  1978.    new.x = new.x * tmp.x - new.y * tmp.y;
  1979.    new.y = new.x * tmp.y - new.y * tmp.x;
  1980.    new.x += floatparm->x;
  1981.    new.y += floatparm->y;
  1982.    return(floatbailout());
  1983. }
  1984. TimsErrorFractal()
  1985. {
  1986. #ifndef XFRACT
  1987.    LCMPLXtrig0(lold,lnew);
  1988.    lnew.x = multiply(lnew.x,ltmp.x,bitshift)-multiply(lnew.y,ltmp.y,bitshift);
  1989.    lnew.y = multiply(lnew.x,ltmp.y,bitshift)-multiply(lnew.y,ltmp.x,bitshift);
  1990.    lnew.x += longparm->x;
  1991.    lnew.y += longparm->y;
  1992.    return(longbailout());
  1993. #endif
  1994. }
  1995.  
  1996. CirclefpFractal()
  1997. {
  1998.    extern int colors;
  1999.    extern int color;
  2000.    int i;
  2001.    i = param[0]*(tempsqrx+tempsqry);
  2002.    color = i&(colors-1);
  2003.    return(1);
  2004. }
  2005. /*
  2006. CirclelongFractal()
  2007. {
  2008.    extern int colors;
  2009.    extern int color;
  2010.    long i;
  2011.    i = multiply(lparm.x,(ltempsqrx+ltempsqry),bitshift);
  2012.    i = i >> bitshift;
  2013.    color = i&(colors-1);
  2014.    return(1);
  2015. }
  2016. */
  2017.  
  2018. /* -------------------------------------------------------------------- */
  2019. /*        Initialization (once per pixel) routines                        */
  2020. /* -------------------------------------------------------------------- */
  2021.  
  2022. #ifdef XFRACT
  2023. /* this code translated to asm - lives in newton.asm */
  2024. /* transform points with reciprocal function */
  2025. void invertz2(_CMPLX *z)
  2026. {
  2027.    z->x = dx0[col]+dx1[row];
  2028.    z->y = dy0[row]+dy1[col];
  2029.    z->x -= f_xcenter; z->y -= f_ycenter;  /* Normalize values to center of circle */
  2030.  
  2031.    tempsqrx = sqr(z->x) + sqr(z->y);  /* Get old radius */
  2032.    if(fabs(tempsqrx) > FLT_MIN)
  2033.       tempsqrx = f_radius / tempsqrx;
  2034.    else
  2035.       tempsqrx = FLT_MAX;   /* a big number, but not TOO big */
  2036.    z->x *= tempsqrx;      z->y *= tempsqrx;     /* Perform inversion */
  2037.    z->x += f_xcenter; z->y += f_ycenter; /* Renormalize */
  2038. }
  2039. #endif
  2040.  
  2041. int long_julia_per_pixel()
  2042. {
  2043. #ifndef XFRACT
  2044.    /* integer julia types */
  2045.    /* lambda */
  2046.    /* barnsleyj1 */
  2047.    /* barnsleyj2 */
  2048.    /* sierpinski */
  2049.    if(invert)
  2050.    {
  2051.       /* invert */
  2052.       invertz2(&old);
  2053.  
  2054.       /* watch out for overflow */
  2055.       if(sqr(old.x)+sqr(old.y) >= 127)
  2056.       {
  2057.      old.x = 8;  /* value to bail out in one iteration */
  2058.      old.y = 8;
  2059.       }
  2060.  
  2061.       /* convert to fudged longs */
  2062.       lold.x = old.x*fudge;
  2063.       lold.y = old.y*fudge;
  2064.    }
  2065.    else
  2066.    {
  2067.       lold.x = lx0[col]+lx1[row];
  2068.       lold.y = ly0[row]+ly1[col];
  2069.    }
  2070.    return(0);
  2071. #else
  2072.    printf("Called long_julia_per_pixel\n");
  2073.    exit(0);
  2074. #endif
  2075. }
  2076.  
  2077. int long_richard8_per_pixel()
  2078. {
  2079. #ifndef XFRACT
  2080.     long_mandel_per_pixel();
  2081.     LCMPLXtrig1(*longparm,ltmp);
  2082.     LCMPLXmult(ltmp,lparm2,ltmp);
  2083.     return(1);
  2084. #endif
  2085. }
  2086.  
  2087. int long_mandel_per_pixel()
  2088. {
  2089. #ifndef XFRACT
  2090.    /* integer mandel types */
  2091.    /* barnsleym1 */
  2092.    /* barnsleym2 */
  2093.    linit.x = lx0[col]+lx1[row];
  2094.  
  2095.    if(invert)
  2096.    {
  2097.       /* invert */
  2098.       invertz2(&init);
  2099.  
  2100.       /* watch out for overflow */
  2101.       if(sqr(init.x)+sqr(init.y) >= 127)
  2102.       {
  2103.      init.x = 8;  /* value to bail out in one iteration */
  2104.      init.y = 8;
  2105.       }
  2106.  
  2107.       /* convert to fudged longs */
  2108.       linit.x = init.x*fudge;
  2109.       linit.y = init.y*fudge;
  2110.    }
  2111.  
  2112.    if(useinitorbit == 1)
  2113.       lold = linitorbit;
  2114.    else
  2115.       lold = linit;
  2116.  
  2117.    lold.x += lparm.x;     /* initial pertubation of parameters set */
  2118.    lold.y += lparm.y;
  2119.    return(1); /* 1st iteration has been done */
  2120. #else
  2121.    printf("Called long_mandel_per_pixel\n");
  2122.    exit(0);
  2123. #endif
  2124. }
  2125.  
  2126. int julia_per_pixel()
  2127. {
  2128.    /* julia */
  2129.  
  2130.    if(invert)
  2131.    {
  2132.       /* invert */
  2133.       invertz2(&old);
  2134.  
  2135.       /* watch out for overflow */
  2136.       if(bitshift <= 24)
  2137.      if (sqr(old.x)+sqr(old.y) >= 127)
  2138.      {
  2139.         old.x = 8;    /* value to bail out in one iteration */
  2140.         old.y = 8;
  2141.      }
  2142.       if(bitshift >  24)
  2143.      if (sqr(old.x)+sqr(old.y) >= 4.0)
  2144.      {
  2145.         old.x = 2;    /* value to bail out in one iteration */
  2146.         old.y = 2;
  2147.      }
  2148.  
  2149.       /* convert to fudged longs */
  2150.       lold.x = old.x*fudge;
  2151.       lold.y = old.y*fudge;
  2152.    }
  2153.    else
  2154.    {
  2155.       lold.x = lx0[col]+lx1[row];
  2156.       lold.y = ly0[row]+ly1[col];
  2157.    }
  2158.  
  2159.    ltempsqrx = multiply(lold.x, lold.x, bitshift);
  2160.    ltempsqry = multiply(lold.y, lold.y, bitshift);
  2161.    ltmp = lold;
  2162.    return(0);
  2163. }
  2164.  
  2165. marks_mandelpwr_per_pixel()
  2166. {
  2167. #ifndef XFRACT
  2168.    mandel_per_pixel();
  2169.    ltmp = lold;
  2170.    ltmp.x -= fudge;
  2171.    LCMPLXpwr(lold,ltmp,ltmp);
  2172.    return(1);
  2173. #endif
  2174. }
  2175.  
  2176. int mandel_per_pixel()
  2177. {
  2178.    /* mandel */
  2179.  
  2180.    if(invert)
  2181.    {
  2182.       invertz2(&init);
  2183.  
  2184.       /* watch out for overflow */
  2185.       if(bitshift <= 24)
  2186.      if (sqr(init.x)+sqr(init.y) >= 127)
  2187.      {
  2188.         init.x = 8;  /* value to bail out in one iteration */
  2189.         init.y = 8;
  2190.      }
  2191.       if(bitshift >  24)
  2192.      if (sqr(init.x)+sqr(init.y) >= 4)
  2193.      {
  2194.         init.x = 2;  /* value to bail out in one iteration */
  2195.         init.y = 2;
  2196.      }
  2197.  
  2198.       /* convert to fudged longs */
  2199.       linit.x = init.x*fudge;
  2200.       linit.y = init.y*fudge;
  2201.    }
  2202.    else
  2203.       linit.x = lx0[col]+lx1[row];
  2204.    switch (fractype)
  2205.      {
  2206.     case MANDELLAMBDA:        /* Critical Value 0.5 + 0.0i  */
  2207.         lold.x = FgHalf;
  2208.         lold.y = 0;
  2209.         break;
  2210.     default:
  2211.         lold = linit;
  2212.         break;
  2213.       }
  2214.  
  2215.    /* alter init value */
  2216.    if(useinitorbit == 1)
  2217.       lold = linitorbit;
  2218.    else if(useinitorbit == 2)
  2219.       lold = linit;
  2220.  
  2221.    if(inside == -60 || inside == -61)
  2222.    {
  2223.       /* kludge to match "Beauty of Fractals" picture since we start
  2224.      Mandelbrot iteration with init rather than 0 */
  2225.       lold.x = lparm.x; /* initial pertubation of parameters set */
  2226.       lold.y = lparm.y;
  2227.       color = -1;
  2228.    }
  2229.    else
  2230.    {
  2231.       lold.x += lparm.x; /* initial pertubation of parameters set */
  2232.       lold.y += lparm.y;
  2233.    }
  2234.    ltmp = linit; /* for spider */
  2235.    ltempsqrx = multiply(lold.x, lold.x, bitshift);
  2236.    ltempsqry = multiply(lold.y, lold.y, bitshift);
  2237.    return(1); /* 1st iteration has been done */
  2238. }
  2239.  
  2240.  
  2241. int marksmandel_per_pixel()
  2242. {
  2243.    /* marksmandel */
  2244.  
  2245.    if(invert)
  2246.    {
  2247.       invertz2(&init);
  2248.  
  2249.       /* watch out for overflow */
  2250.       if(sqr(init.x)+sqr(init.y) >= 127)
  2251.       {
  2252.      init.x = 8;  /* value to bail out in one iteration */
  2253.      init.y = 8;
  2254.       }
  2255.  
  2256.       /* convert to fudged longs */
  2257.       linit.x = init.x*fudge;
  2258.       linit.y = init.y*fudge;
  2259.    }
  2260.    else
  2261.       linit.x = lx0[col]+lx1[row];
  2262.  
  2263.    if(useinitorbit == 1)
  2264.       lold = linitorbit;
  2265.    else
  2266.       lold = linit;
  2267.  
  2268.    lold.x += lparm.x;     /* initial pertubation of parameters set */
  2269.    lold.y += lparm.y;
  2270.  
  2271.    if(c_exp > 3)
  2272.       lcpower(&lold,c_exp-1,&lcoefficient,bitshift);
  2273.    else if(c_exp == 3) {
  2274.       lcoefficient.x = multiply(lold.x, lold.x, bitshift)
  2275.      - multiply(lold.y, lold.y, bitshift);
  2276.       lcoefficient.y = multiply(lold.x, lold.y, bitshiftless1);
  2277.    }
  2278.    else if(c_exp == 2)
  2279.       lcoefficient = lold;
  2280.    else if(c_exp < 2) {
  2281.       lcoefficient.x = 1L << bitshift;
  2282.       lcoefficient.y = 0L;
  2283.    }
  2284.  
  2285.    ltempsqrx = multiply(lold.x, lold.x, bitshift);
  2286.    ltempsqry = multiply(lold.y, lold.y, bitshift);
  2287.    return(1); /* 1st iteration has been done */
  2288. }
  2289.  
  2290. int marksmandelfp_per_pixel()
  2291. {
  2292.    /* marksmandel */
  2293.  
  2294.    if(invert)
  2295.       invertz2(&init);
  2296.    else
  2297.       init.x = dx0[col]+dx1[row];
  2298.  
  2299.    if(useinitorbit == 1)
  2300.       old = initorbit;
  2301.    else
  2302.       old = init;
  2303.  
  2304.    old.x += parm.x;     /* initial pertubation of parameters set */
  2305.    old.y += parm.y;
  2306.  
  2307.    tempsqrx = sqr(old.x);
  2308.    tempsqry = sqr(old.y);
  2309.  
  2310.    if(c_exp > 3)
  2311.       cpower(&old,c_exp-1,&coefficient);
  2312.    else if(c_exp == 3) {
  2313.       coefficient.x = tempsqrx - tempsqry;
  2314.       coefficient.y = old.x * old.y * 2;
  2315.    }
  2316.    else if(c_exp == 2)
  2317.       coefficient = old;
  2318.    else if(c_exp < 2) {
  2319.       coefficient.x = 1.0;
  2320.       coefficient.y = 0.0;
  2321.    }
  2322.  
  2323.    return(1); /* 1st iteration has been done */
  2324. }
  2325.  
  2326. marks_mandelpwrfp_per_pixel()
  2327. {
  2328.    mandelfp_per_pixel();
  2329.    tmp = old;
  2330.    tmp.x -= 1;
  2331.    CMPLXpwr(old,tmp,tmp);
  2332.    return(1);
  2333. }
  2334.  
  2335. int mandelfp_per_pixel()
  2336. {
  2337.    /* floating point mandelbrot */
  2338.    /* mandelfp */
  2339.  
  2340.    if(invert)
  2341.       invertz2(&init);
  2342.    else
  2343.       init.x = dx0[col]+dx1[row];
  2344.     switch (fractype)
  2345.       {
  2346.     case MAGNET2M:
  2347.         FloatPreCalcMagnet2();
  2348.     case MAGNET1M:         /* Crit Val Zero both, but neither   */
  2349.         old.x = old.y = 0.0; /* is of the form f(Z,C) = Z*g(Z)+C  */
  2350.         break;
  2351.     case MANDELLAMBDAFP:        /* Critical Value 0.5 + 0.0i  */
  2352.         old.x = 0.5;
  2353.         old.y = 0.0;
  2354.         break;
  2355.     default:
  2356.         old = init;
  2357.         break;
  2358.       }
  2359.  
  2360.    /* alter init value */
  2361.    if(useinitorbit == 1)
  2362.       old = initorbit;
  2363.    else if(useinitorbit == 2)
  2364.       old = init;
  2365.  
  2366.    if(inside == -60 || inside == -61)
  2367.    {
  2368.       /* kludge to match "Beauty of Fractals" picture since we start
  2369.      Mandelbrot iteration with init rather than 0 */
  2370.       old.x = parm.x; /* initial pertubation of parameters set */
  2371.       old.y = parm.y;
  2372.       color = -1;
  2373.    }
  2374.    else
  2375.    {
  2376.      old.x += parm.x;
  2377.      old.y += parm.y;
  2378.    }
  2379.    tmp = init; /* for spider */
  2380.    tempsqrx = sqr(old.x);  /* precalculated value for regular Mandelbrot */
  2381.    tempsqry = sqr(old.y);
  2382.    return(1); /* 1st iteration has been done */
  2383. }
  2384.  
  2385. int juliafp_per_pixel()
  2386. {
  2387.    /* floating point julia */
  2388.    /* juliafp */
  2389.    if(invert)
  2390.       invertz2(&old);
  2391.    else
  2392.    {
  2393.      old.x = dx0[col]+dx1[row];
  2394.      old.y = dy0[row]+dy1[col];
  2395.    }
  2396.    tempsqrx = sqr(old.x);  /* precalculated value for regular Julia */
  2397.    tempsqry = sqr(old.y);
  2398.    tmp = old;
  2399.    return(0);
  2400. }
  2401.  
  2402. int MPCjulia_per_pixel()
  2403. {
  2404. #ifndef XFRACT
  2405.    /* floating point julia */
  2406.    /* juliafp */
  2407.    if(invert)
  2408.       invertz2(&old);
  2409.    else
  2410.    {
  2411.      old.x = dx0[col]+dx1[row];
  2412.      old.y = dy0[row]+dy1[col];
  2413.    }
  2414.    mpcold.x = *pd2MP(old.x);
  2415.    mpcold.y = *pd2MP(old.y);
  2416.    return(0);
  2417. #endif
  2418. }
  2419.  
  2420. otherrichard8fp_per_pixel()
  2421. {
  2422.     othermandelfp_per_pixel();
  2423.     CMPLXtrig1(*floatparm,tmp);
  2424.     CMPLXmult(tmp,parm2,tmp);
  2425.     return(1);
  2426. }
  2427.  
  2428. int othermandelfp_per_pixel()
  2429. {
  2430.    if(invert)
  2431.       invertz2(&init);
  2432.    else
  2433.       init.x = dx0[col]+dx1[row];
  2434.  
  2435.    if(useinitorbit == 1)
  2436.       old = initorbit;
  2437.    else
  2438.       old = init;
  2439.  
  2440.    old.x += parm.x;     /* initial pertubation of parameters set */
  2441.    old.y += parm.y;
  2442.  
  2443.    return(1); /* 1st iteration has been done */
  2444. }
  2445.  
  2446. int MPCHalley_per_pixel()
  2447. {
  2448. #ifndef XFRACT
  2449.    /* MPC halley */
  2450.    if(invert)
  2451.       invertz2(&init);
  2452.    else
  2453.      init.x = dx0[col]+dx1[row];
  2454.  
  2455.    mpcold.x = *pd2MP(init.x);
  2456.    mpcold.y = *pd2MP(init.y);
  2457.  
  2458.    return(0);
  2459. #endif
  2460. }
  2461.  
  2462. int Halley_per_pixel()
  2463. {
  2464.    if(invert)
  2465.       invertz2(&init);
  2466.    else
  2467.       init.x = dx0[col]+dx1[row];
  2468.  
  2469.    old = init;
  2470.  
  2471.    return(0); /* 1st iteration is not done */
  2472. }
  2473.  
  2474. int otherjuliafp_per_pixel()
  2475. {
  2476.    if(invert)
  2477.       invertz2(&old);
  2478.    else
  2479.    {
  2480.       old.x = dx0[col]+dx1[row];
  2481.       old.y = dy0[row]+dy1[col];
  2482.    }
  2483.    return(0);
  2484. }
  2485.  
  2486. #if 0
  2487. #define Q0 .113
  2488. #define Q1 .01
  2489. #else
  2490. #define Q0 0
  2491. #define Q1 0
  2492. #endif
  2493.  
  2494. int quaternionjulfp_per_pixel()
  2495. {
  2496.    old.x = dx0[col]+dx1[row];
  2497.    old.y = dy0[row]+dy1[col];
  2498.    floatparm->x = param[4];
  2499.    floatparm->y = param[5];
  2500.    qc  = param[0];
  2501.    qci = param[1];
  2502.    qcj = param[2];
  2503.    qck = param[3];
  2504.    return(0);
  2505. }
  2506.  
  2507. int quaternionfp_per_pixel()
  2508. {
  2509.    old.x = 0;
  2510.    old.y = 0;
  2511.    floatparm->x = 0;
  2512.    floatparm->y = 0;
  2513.    qc  = dx0[col]+dx1[row];
  2514.    qci = dy0[row]+dy1[col];
  2515.    qcj = param[2];
  2516.    qck = param[3];
  2517.    return(0);
  2518. }
  2519.  
  2520. int trigmandelfp_per_pixel()
  2521. {
  2522.    if(invert)
  2523.       invertz2(&init);
  2524.    else
  2525.       init.x = dx0[col]+dx1[row];
  2526.  
  2527.    if(useinitorbit == 1)
  2528.       old = initorbit;
  2529.    else
  2530.       old = init;
  2531.  
  2532.    old.x += parm.x;     /* initial pertubation of parameters set */
  2533.    old.y += parm.y;
  2534.    CMPLXtrig0(old,old);
  2535.    return(1); /* 1st iteration has been done */
  2536. }
  2537.  
  2538. int trigjuliafp_per_pixel()
  2539. {
  2540.    /* for tetrated types */
  2541.    if(invert)
  2542.       invertz2(&old);
  2543.    else
  2544.    {
  2545.       old.x = dx0[col]+dx1[row];
  2546.       old.y = dy0[row]+dy1[col];
  2547.    }
  2548.    CMPLXtrig0(old,old);
  2549.    return(0);
  2550. }
  2551.  
  2552. int trigXtrigmandelfp_per_pixel()
  2553. {
  2554.    if(invert)
  2555.       invertz2(&init);
  2556.    else
  2557.       init.x = dx0[col]+dx1[row];
  2558.  
  2559.    if(useinitorbit == 1)
  2560.       old = initorbit;
  2561.    else
  2562.       old = init;
  2563.  
  2564.    old.x += parm.x;     /* initial pertubation of parameters set */
  2565.    old.y += parm.y;
  2566.    CMPLXtrig0(old,tmp);
  2567.    CMPLXtrig1(old,tmp2);
  2568.    CMPLXmult(tmp,tmp2,old);
  2569.    return(1); /* 1st iteration has been done */
  2570. }
  2571.  
  2572. int trigXtrigjuliafp_per_pixel()
  2573. {
  2574.    if(invert)
  2575.       invertz2(&old);
  2576.    else
  2577.    {
  2578.       old.x = dx0[col]+dx1[row];
  2579.       old.y = dy0[row]+dy1[col];
  2580.    }
  2581.    CMPLXtrig0(old,tmp);
  2582.    CMPLXtrig1(old,tmp2);
  2583.    CMPLXmult(tmp,tmp2,old);
  2584.    return(0);
  2585. }
  2586.  
  2587. int MarksCplxMandperp(void)
  2588. {
  2589.    if(invert)
  2590.       invertz2(&init);
  2591.    else
  2592.       init.x = dx0[col]+dx1[row];
  2593.    old.x = init.x + parm.x; /* initial pertubation of parameters set */
  2594.    old.y = init.y + parm.y;
  2595.    tempsqrx = sqr(old.x);  /* precalculated value */
  2596.    tempsqry = sqr(old.y);
  2597.    Coefficient = ComplexPower(init, pwr);
  2598.    return(1);
  2599. }
  2600.  
  2601. int long_phoenix_per_pixel()
  2602. {
  2603. #ifndef XFRACT
  2604.    if(invert)
  2605.    {
  2606.       /* invert */
  2607.       invertz2(&old);
  2608.  
  2609.       /* watch out for overflow */
  2610.       if(sqr(old.x)+sqr(old.y) >= 127)
  2611.       {
  2612.      old.x = 8;  /* value to bail out in one iteration */
  2613.      old.y = 8;
  2614.       }
  2615.  
  2616.       /* convert to fudged longs */
  2617.       lold.x = old.x*fudge;
  2618.       lold.y = old.y*fudge;
  2619.    }
  2620.    else
  2621.    {
  2622.       lold.x = lx0[col]+lx1[row];
  2623.       lold.y = ly0[row]+ly1[col];
  2624.    }
  2625.    ltempsqrx = multiply(lold.x, lold.x, bitshift);
  2626.    ltempsqry = multiply(lold.y, lold.y, bitshift);
  2627.    ltmp2.x = 0; /* use ltmp2 as the complex Y value */
  2628.    ltmp2.y = 0;
  2629.    return(0);
  2630. #else
  2631.    printf("Called long_phoenix_per_pixel\n");
  2632.    exit(0);
  2633. #endif
  2634. }
  2635.  
  2636. int phoenix_per_pixel()
  2637. {
  2638.    if(invert)
  2639.       invertz2(&old);
  2640.    else
  2641.    {
  2642.       old.x = dx0[col]+dx1[row];
  2643.       old.y = dy0[row]+dy1[col];
  2644.    }
  2645.    tempsqrx = sqr(old.x);  /* precalculated value */
  2646.    tempsqry = sqr(old.y);
  2647.    tmp2.x = 0; /* use tmp2 as the complex Y value */
  2648.    tmp2.y = 0;
  2649.    return(0);
  2650. }
  2651. int long_mandphoenix_per_pixel()
  2652. {
  2653. #ifndef XFRACT
  2654.    linit.x = lx0[col]+lx1[row];
  2655.  
  2656.    if(invert)
  2657.    {
  2658.       /* invert */
  2659.       invertz2(&init);
  2660.  
  2661.       /* watch out for overflow */
  2662.       if(sqr(init.x)+sqr(init.y) >= 127)
  2663.       {
  2664.      init.x = 8;  /* value to bail out in one iteration */
  2665.      init.y = 8;
  2666.       }
  2667.  
  2668.       /* convert to fudged longs */
  2669.       linit.x = init.x*fudge;
  2670.       linit.y = init.y*fudge;
  2671.    }
  2672.  
  2673.    if(useinitorbit == 1)
  2674.       lold = linitorbit;
  2675.    else
  2676.       lold = linit;
  2677.  
  2678.    lold.x += lparm.x;     /* initial pertubation of parameters set */
  2679.    lold.y += lparm.y;
  2680.    ltempsqrx = multiply(lold.x, lold.x, bitshift);
  2681.    ltempsqry = multiply(lold.y, lold.y, bitshift);
  2682.    ltmp2.x = 0;
  2683.    ltmp2.y = 0;
  2684.    return(1); /* 1st iteration has been done */
  2685. #else
  2686.    printf("Called long_mandphoenix_per_pixel\n");
  2687.    exit(0);
  2688. #endif
  2689. }
  2690. int mandphoenix_per_pixel()
  2691. {
  2692.    if(invert)
  2693.       invertz2(&init);
  2694.    else
  2695.       init.x = dx0[col]+dx1[row];
  2696.  
  2697.    if(useinitorbit == 1)
  2698.       old = initorbit;
  2699.    else
  2700.       old = init;
  2701.  
  2702.    old.x += parm.x;     /* initial pertubation of parameters set */
  2703.    old.y += parm.y;
  2704.    tempsqrx = sqr(old.x);  /* precalculated value */
  2705.    tempsqry = sqr(old.y);
  2706.    tmp2.x = 0;
  2707.    tmp2.y = 0;
  2708.    return(1); /* 1st iteration has been done */
  2709. }
  2710.  
  2711. QuaternionFPFractal()
  2712. {
  2713.    double a0,a1,a2,a3,n0,n1,n2,n3;
  2714.    a0 = old.x;
  2715.    a1 = old.y;
  2716.    a2 = floatparm->x;
  2717.    a3 = floatparm->y;
  2718.  
  2719.    n0 = a0*a0-a1*a1-a2*a2-a3*a3 + qc;
  2720.    n1 = 2*a0*a1 + qci;
  2721.    n2 = 2*a0*a2 + qcj;
  2722.    n3 = 2*a0*a3 + qck;
  2723.    /* Check bailout */
  2724.    magnitude = a0*a0+a1*a1+a2*a2+a3*a3;
  2725.    if (magnitude>rqlim) {
  2726.        return 1;
  2727.    }
  2728.    old.x = n0;
  2729.    old.y = n1;
  2730.    floatparm->x = n2;
  2731.    floatparm->y = n3;
  2732.    return(0);
  2733. }
  2734.  
  2735. HyperComplex1FPFractal()
  2736. {
  2737.    double n0,n1,n2,n3;
  2738.  
  2739.    n0 = sqr(old.x)-sqr(old.y)-sqr(floatparm->x)+sqr(floatparm->y) + qc;
  2740.    n1 = 2*old.x*old.y - 2*floatparm->x*floatparm->y + qci;
  2741.    n2 = 2*old.x*floatparm->x - 2*old.y*floatparm->y + qcj;
  2742.    n3 = 2*old.x*floatparm->y + 2*old.y*floatparm->x + qck;
  2743.    
  2744.    old.x = n0;
  2745.    old.y = n1;
  2746.    floatparm->x = n2;
  2747.    floatparm->y = n3;
  2748.  
  2749.    /* Check bailout */
  2750.    magnitude = sqr(old.x)+sqr(old.y)+sqr(floatparm->x)+sqr(floatparm->y);
  2751.    if (magnitude>rqlim) {
  2752.        return 1;
  2753.    }
  2754.    return(0);
  2755. }
  2756.  
  2757. HyperComplexFPFractal()
  2758. {
  2759.    double n0,n1,n2,n3;
  2760.    _HCMPLX hold, hnew;
  2761.    hold.x = old.x;
  2762.    hold.y = old.y;
  2763.    hold.z = floatparm->x;
  2764.    hold.t = floatparm->y;
  2765.  
  2766. /*   HComplexSqr(&hold,&hnew); */
  2767.    HComplexTrig0(&hold,&hnew);
  2768.    
  2769.    hnew.x += qc;
  2770.    hnew.y += qci;
  2771.    hnew.z += qcj;
  2772.    hnew.t += qck;
  2773.    
  2774.    old.x = hnew.x;
  2775.    old.y = hnew.y;
  2776.    floatparm->x = hnew.z;
  2777.    floatparm->y = hnew.t;
  2778.  
  2779.    /* Check bailout */
  2780.    magnitude = sqr(old.x)+sqr(old.y)+sqr(floatparm->x)+sqr(floatparm->y);
  2781.    if (magnitude>rqlim) {
  2782.        return 1;
  2783.    }
  2784.    return(0);
  2785. }
  2786.  
  2787. /* -------------------------------------------------------------------- */
  2788. /*        Setup (once per fractal image) routines         */
  2789. /* -------------------------------------------------------------------- */
  2790.  
  2791. MandelSetup()        /* Mandelbrot Routine */
  2792. {
  2793.    if (debugflag != 90 && ! invert && decomp[0] == 0 && rqlim <= 4.0
  2794.        && bitshift == 29 && potflag == 0
  2795.        && biomorph == -1 && inside > -59 && outside >= -1
  2796.        && useinitorbit != 1 && using_jiim == 0)
  2797.       calctype = calcmand; /* the normal case - use CALCMAND */
  2798.    else
  2799.    {
  2800.       /* special case: use the main processing loop */
  2801.       calctype = StandardFractal;
  2802.       longparm = &linit;
  2803.    }
  2804.    return(1);
  2805. }
  2806.  
  2807. JuliaSetup()        /* Julia Routine */
  2808. {
  2809.    if (debugflag != 90 && ! invert && decomp[0] == 0 && rqlim <= 4.0
  2810.        && bitshift == 29 && potflag == 0
  2811.        && biomorph == -1 && inside > -59 && outside >= -1
  2812.        && !finattract && using_jiim == 0)
  2813.       calctype = calcmand; /* the normal case - use CALCMAND */
  2814.    else
  2815.    {
  2816.       /* special case: use the main processing loop */
  2817.       calctype = StandardFractal;
  2818.       longparm = &lparm;
  2819.       get_julia_attractor (0.0, 0.0);    /* another attractor? */
  2820.    }
  2821.    return(1);
  2822. }
  2823.  
  2824. NewtonSetup()        /* Newton/NewtBasin Routines */
  2825. {
  2826.    int i;
  2827.    extern int basin;
  2828.    extern int fpu;
  2829.    if (debugflag != 1010)
  2830.    {
  2831.       if(fpu != 0)
  2832.       {
  2833.      if(fractype == MPNEWTON)
  2834.         fractype = NEWTON;
  2835.      else if(fractype == MPNEWTBASIN)
  2836.         fractype = NEWTBASIN;
  2837.       }
  2838.       else
  2839.       {
  2840.      if(fractype == NEWTON)
  2841.            fractype = MPNEWTON;
  2842.      else if(fractype == NEWTBASIN)
  2843.            fractype = MPNEWTBASIN;
  2844.       }
  2845.       curfractalspecific = &fractalspecific[fractype];
  2846.    }
  2847.  
  2848.    /* set up table of roots of 1 along unit circle */
  2849.    degree = (int)parm.x;
  2850.    if(degree < 2)
  2851.       degree = 3;   /* defaults to 3, but 2 is possible */
  2852.    root = 1;
  2853.  
  2854.    /* precalculated values */
  2855.    roverd    = (double)root / (double)degree;
  2856.    d1overd    = (double)(degree - 1) / (double)degree;
  2857.    maxcolor    = 0;
  2858.    threshold    = .3*PI/degree; /* less than half distance between roots */
  2859. #ifndef XFRACT
  2860.    if (fractype == MPNEWTON || fractype == MPNEWTBASIN) {
  2861.       mproverd       = *pd2MP(roverd);
  2862.       mpd1overd    = *pd2MP(d1overd);
  2863.       mpthreshold  = *pd2MP(threshold);
  2864.       mpone       = *pd2MP(1.0);
  2865.    }
  2866. #endif
  2867.  
  2868.    floatmin = FLT_MIN;
  2869.    floatmax = FLT_MAX;
  2870.  
  2871.    basin = 0;
  2872.    if(roots != staticroots) {
  2873.       free(roots);
  2874.       roots = staticroots;
  2875.    }
  2876.  
  2877.    if (fractype==NEWTBASIN)
  2878.    {
  2879.       if(parm.y)
  2880.      basin = 2; /*stripes */
  2881.       else
  2882.      basin = 1;
  2883.       if(degree > 16)
  2884.       {
  2885.      if((roots=(_CMPLX *)malloc(degree*sizeof(_CMPLX)))==NULL)
  2886.      {
  2887.         roots = staticroots;
  2888.         degree = 16;
  2889.      }
  2890.       }
  2891.       else
  2892.      roots = staticroots;
  2893.  
  2894.       /* list of roots to discover where we converged for newtbasin */
  2895.       for(i=0;i<degree;i++)
  2896.       {
  2897.      roots[i].x = cos(i*twopi/(double)degree);
  2898.      roots[i].y = sin(i*twopi/(double)degree);
  2899.       }
  2900.    }
  2901. #ifndef XFRACT
  2902.    else if (fractype==MPNEWTBASIN)
  2903.    {
  2904.      if(parm.y)
  2905.      basin = 2; /*stripes */
  2906.       else
  2907.      basin = 1;
  2908.  
  2909.       if(degree > 16)
  2910.       {
  2911.      if((MPCroots=(struct MPC *)malloc(degree*sizeof(struct MPC)))==NULL)
  2912.      {
  2913.         MPCroots = (struct MPC *)staticroots;
  2914.         degree = 16;
  2915.      }
  2916.       }
  2917.       else
  2918.      MPCroots = (struct MPC *)staticroots;
  2919.  
  2920.       /* list of roots to discover where we converged for newtbasin */
  2921.       for(i=0;i<degree;i++)
  2922.       {
  2923.      MPCroots[i].x = *pd2MP(cos(i*twopi/(double)degree));
  2924.      MPCroots[i].y = *pd2MP(sin(i*twopi/(double)degree));
  2925.       }
  2926.    }
  2927. #endif
  2928.  
  2929.    param[0] = (double)degree; /* JCO 7/1/92 */
  2930.    if (degree%4 == 0)
  2931.       symmetry = XYAXIS;
  2932.    else
  2933.       symmetry = XAXIS;
  2934.  
  2935.    calctype=StandardFractal;
  2936. #ifndef XFRACT
  2937.    if (fractype == MPNEWTON || fractype == MPNEWTBASIN)
  2938.       setMPfunctions();
  2939. #endif
  2940.    return(1);
  2941. }
  2942.  
  2943.  
  2944. StandaloneSetup()
  2945. {
  2946.    timer(0,curfractalspecific->calctype);
  2947.    return(0);        /* effectively disable solid-guessing */
  2948. }
  2949.  
  2950. UnitySetup()
  2951. {
  2952.    periodicitycheck = 0;
  2953.    FgOne = (1L << bitshift);
  2954.    FgTwo = FgOne + FgOne;
  2955.    return(1);
  2956. }
  2957.  
  2958. MandelfpSetup()
  2959. {
  2960.    c_exp = param[2];
  2961.    pwr.x = param[2] - 1.0;
  2962.    pwr.y = param[3];
  2963.    floatparm = &init;
  2964.    switch (fractype)
  2965.    {
  2966.    case MARKSMANDELFP:
  2967.       if(c_exp < 1)
  2968.          c_exp = 1;
  2969.       if(!(c_exp & 1))
  2970.          symmetry = XYAXIS_NOPARM;    /* odd exponents */
  2971.       if(c_exp & 1)
  2972.          symmetry = XAXIS_NOPARM;
  2973.       break;
  2974.    case MANDELFP:
  2975.     /*
  2976.        floating point code could probably be altered to handle many of
  2977.        the situations that otherwise are using StandardFractal().
  2978.        calcmandfp() can currently handle invert, any rqlim, potflag
  2979.        zmag, epsilon cross, and all the current outside options
  2980.                              Wes Loewer 11/03/91
  2981.     */
  2982.     if (debugflag != 90
  2983.         && !distest
  2984.         && decomp[0] == 0
  2985.         && biomorph == -1
  2986.         && (inside >= -1 || inside == -59 || inside == -100)
  2987.         /* uncomment this next line if more outside options are added */
  2988.         /* && outside >= -5 */
  2989.         && useinitorbit != 1
  2990.         && using_jiim == 0)
  2991.     {
  2992.        calctype = calcmandfp; /* the normal case - use calcmandfp */
  2993.        calcmandfpasmstart();
  2994.     }
  2995.     else
  2996.     {
  2997.        /* special case: use the main processing loop */
  2998.        calctype = StandardFractal;
  2999.     }
  3000.     break;
  3001.    case FPMANDELZPOWER:
  3002.       if((double)c_exp == param[2] && c_exp & 1) /* odd exponents */
  3003.          symmetry = XYAXIS_NOPARM;
  3004.       if(param[3] != 0)
  3005.          symmetry = NOSYM;
  3006.       if(param[3] == 0.0 && debugflag != 6000 && (double)c_exp == param[2])
  3007.           fractalspecific[fractype].orbitcalc = floatZpowerFractal;
  3008.       else
  3009.           fractalspecific[fractype].orbitcalc = floatCmplxZpowerFractal;
  3010.       break;
  3011.    case MAGNET1M:
  3012.    case MAGNET2M:
  3013.       attr[0].x = 1.0;        /* 1.0 + 0.0i always attracts */
  3014.       attr[0].y = 0.0;        /* - both MAGNET1 and MAGNET2 */
  3015.       attrperiod[0] = 1;
  3016.       attractors = 1;
  3017.       break;
  3018.    case SPIDERFP:
  3019.       if(periodicitycheck==1) /* if not user set */
  3020.      periodicitycheck=4;
  3021.       break;
  3022.    case MANDELEXP:
  3023.       symmetry = XAXIS_NOPARM;
  3024.       break;
  3025. /* Added to account for symmetry in manfn+exp and manfn+zsqrd */
  3026. /*     JCO 2/29/92 */
  3027.    case FPMANTRIGPLUSEXP:
  3028.    case FPMANTRIGPLUSZSQRD:
  3029.      if(parm.y == 0.0)
  3030.         symmetry = XAXIS;
  3031.      else
  3032.         symmetry = NOSYM;
  3033.      if ((trigndx[0] == LOG) || (trigndx[0] == 14)) /* LOG or FLIP */
  3034.         symmetry = NOSYM;
  3035.       break;
  3036.    case QUATFP:
  3037.    case HYPERCMPLXFP:
  3038.       floatparm = &tmp;
  3039.       attractors = 0;
  3040.       periodicitycheck = 0;
  3041.       break;
  3042.    default:
  3043.       break;
  3044.    }
  3045.    return(1);
  3046. }
  3047.  
  3048. JuliafpSetup()
  3049. {
  3050.    c_exp = param[2];
  3051.    floatparm = &parm;
  3052.    if(fractype==COMPLEXMARKSJUL)
  3053.    {
  3054.       pwr.x = param[2] - 1.0;
  3055.       pwr.y = param[3];
  3056.       Coefficient = ComplexPower(*floatparm, pwr);
  3057.    }
  3058.    switch (fractype)
  3059.    {
  3060.    case JULIAFP:
  3061.     /*
  3062.        floating point code could probably be altered to handle many of
  3063.        the situations that otherwise are using StandardFractal().
  3064.        calcmandfp() can currently handle invert, any rqlim, potflag
  3065.        zmag, epsilon cross, and all the current outside options
  3066.                              Wes Loewer 11/03/91
  3067.     */
  3068.     if (debugflag != 90
  3069.         && !distest
  3070.         && decomp[0] == 0
  3071.         && biomorph == -1
  3072.         && (inside >= -1 || inside == -59 || inside == -100)
  3073.         /* uncomment this next line if more outside options are added */
  3074.         /* && outside >= -5 */
  3075.         && useinitorbit != 1
  3076.         && !finattract
  3077.         && using_jiim == 0)
  3078.     {
  3079.        calctype = calcmandfp; /* the normal case - use calcmandfp */
  3080.        calcmandfpasmstart();
  3081.     }
  3082.     else
  3083.     {
  3084.        /* special case: use the main processing loop */
  3085.        calctype = StandardFractal;
  3086.        get_julia_attractor (0.0, 0.0);   /* another attractor? */
  3087.     }
  3088.     break;
  3089.    case FPJULIAZPOWER:
  3090.       if((c_exp & 1) || param[3] != 0.0 || (double)c_exp != param[2] )
  3091.          symmetry = NOSYM;
  3092.       if(param[3] == 0.0 && debugflag != 6000 && (double)c_exp == param[2])
  3093.           fractalspecific[fractype].orbitcalc = floatZpowerFractal;
  3094.       else
  3095.           fractalspecific[fractype].orbitcalc = floatCmplxZpowerFractal;
  3096.       get_julia_attractor (param[0], param[1]);    /* another attractor? */
  3097.       break;
  3098.    case MAGNET2J:
  3099.       FloatPreCalcMagnet2();
  3100.    case MAGNET1J:
  3101.       attr[0].x = 1.0;        /* 1.0 + 0.0i always attracts */
  3102.       attr[0].y = 0.0;        /* - both MAGNET1 and MAGNET2 */
  3103.       attrperiod[0] = 1;
  3104.       attractors = 1;
  3105.       get_julia_attractor (0.0, 0.0);    /* another attractor? */
  3106.       break;
  3107.    case LAMBDAFP:
  3108.       get_julia_attractor (0.0, 0.0);    /* another attractor? */
  3109.       get_julia_attractor (0.5, 0.0);    /* another attractor? */
  3110.       break;
  3111.    case LAMBDAEXP:
  3112.       if(parm.y == 0.0)
  3113.      symmetry=XAXIS;
  3114.       get_julia_attractor (0.0, 0.0);    /* another attractor? */
  3115.       break;
  3116. /* Added to account for symmetry in julfn+exp and julfn+zsqrd */
  3117. /*     JCO 2/29/92 */
  3118.    case FPJULTRIGPLUSEXP:
  3119.    case FPJULTRIGPLUSZSQRD:
  3120.      if(parm.y == 0.0)
  3121.         symmetry = XAXIS;
  3122.      else
  3123.         symmetry = NOSYM;
  3124.      if ((trigndx[0] == LOG) || (trigndx[0] == 14)) /* LOG or FLIP */
  3125.         symmetry = NOSYM;
  3126.       get_julia_attractor (0.0, 0.0);    /* another attractor? */
  3127.       break;
  3128.    case HYPERCMPLXJFP:
  3129.       if(trigndx[0] != SQR)
  3130.          symmetry=NOSYM;
  3131.    case QUATJULFP:
  3132.       attractors = 0;    /* attractors broken since code checks r,i not j,k */
  3133.       periodicitycheck = 0;
  3134.       if(param[4] != 0.0 || param[5] != 0)
  3135.          symmetry = NOSYM;
  3136.       break;
  3137.    default:
  3138.       get_julia_attractor (0.0, 0.0);    /* another attractor? */
  3139.       break;
  3140.    }
  3141.    return(1);
  3142. }
  3143.  
  3144. MandellongSetup()
  3145. {
  3146.    FgHalf = fudge >> 1;
  3147.    c_exp = param[2];
  3148.    if(fractype==MARKSMANDEL && c_exp < 1)
  3149.       c_exp = 1;
  3150.    if(fractype==LMANDELZPOWER && c_exp < 1)
  3151.       c_exp = 1;
  3152.    if((fractype==MARKSMANDEL   && !(c_exp & 1)) ||
  3153.        (fractype==LMANDELZPOWER && c_exp & 1))
  3154.       symmetry = XYAXIS_NOPARM;    /* odd exponents */
  3155.    if((fractype==MARKSMANDEL && (c_exp & 1)) || fractype==LMANDELEXP)
  3156.       symmetry = XAXIS_NOPARM;
  3157.    if(fractype==SPIDER && periodicitycheck==1)
  3158.       periodicitycheck=4;
  3159.    longparm = &linit;
  3160.    if(fractype==LMANDELZPOWER)
  3161.    {
  3162.       if(param[3] == 0.0 && debugflag != 6000  && (double)c_exp == param[2])
  3163.           fractalspecific[fractype].orbitcalc = longZpowerFractal;
  3164.       else
  3165.           fractalspecific[fractype].orbitcalc = longCmplxZpowerFractal;
  3166.       if(param[3] != 0 || (double)c_exp != param[2] )
  3167.          symmetry = NOSYM;
  3168.     }
  3169. /* Added to account for symmetry in manfn+exp and manfn+zsqrd */
  3170. /*     JCO 2/29/92 */
  3171.    if((fractype==LMANTRIGPLUSEXP)||(fractype==LMANTRIGPLUSZSQRD))
  3172.    {
  3173.      if(parm.y == 0.0)
  3174.         symmetry = XAXIS;
  3175.      else
  3176.         symmetry = NOSYM;
  3177.      if ((trigndx[0] == LOG) || (trigndx[0] == 14)) /* LOG or FLIP */
  3178.         symmetry = NOSYM;
  3179.    }
  3180.    return(1);
  3181. }
  3182.  
  3183. JulialongSetup()
  3184. {
  3185.    c_exp = param[2];
  3186.    longparm = &lparm;
  3187.    switch (fractype)
  3188.    {
  3189.    case LJULIAZPOWER:
  3190.       if((c_exp & 1) || param[3] != 0.0 || (double)c_exp != param[2])
  3191.          symmetry = NOSYM;
  3192.       else if(c_exp < 1)
  3193.          c_exp = 1;
  3194.       if(param[3] == 0.0 && debugflag != 6000 && (double)c_exp == param[2])
  3195.           fractalspecific[fractype].orbitcalc = longZpowerFractal;
  3196.       else
  3197.           fractalspecific[fractype].orbitcalc = longCmplxZpowerFractal;
  3198.       break;
  3199.    case LAMBDA:
  3200.       get_julia_attractor (0.0, 0.0);    /* another attractor? */
  3201.       get_julia_attractor (0.5, 0.0);    /* another attractor? */
  3202.       break;
  3203.    case LLAMBDAEXP:
  3204.       if(lparm.y == 0)
  3205.      symmetry = XAXIS;
  3206.       break;
  3207. /* Added to account for symmetry in julfn+exp and julfn+zsqrd */
  3208. /*     JCO 2/29/92 */
  3209.    case LJULTRIGPLUSEXP:
  3210.    case LJULTRIGPLUSZSQRD:
  3211.      if(parm.y == 0.0)
  3212.         symmetry = XAXIS;
  3213.      else
  3214.         symmetry = NOSYM;
  3215.      if ((trigndx[0] == LOG) || (trigndx[0] == 14)) /* LOG or FLIP */
  3216.         symmetry = NOSYM;
  3217.       get_julia_attractor (0.0, 0.0);    /* another attractor? */
  3218.       break;
  3219.    default:
  3220.       get_julia_attractor (0.0, 0.0);    /* another attractor? */
  3221.       break;
  3222.    }
  3223.    return(1);
  3224. }
  3225.  
  3226. TrigPlusSqrlongSetup()
  3227. {
  3228.    curfractalspecific->per_pixel =  julia_per_pixel;
  3229.    curfractalspecific->orbitcalc =  TrigPlusSqrFractal;
  3230.    if(lparm.x == fudge && lparm.y == 0L && lparm2.y == 0L && debugflag != 90)
  3231.    {
  3232.       if(lparm2.x == fudge)       /* Scott variant */
  3233.      curfractalspecific->orbitcalc =  ScottTrigPlusSqrFractal;
  3234.       else if(lparm2.x == -fudge)  /* Skinner variant */
  3235.      curfractalspecific->orbitcalc =  SkinnerTrigSubSqrFractal;
  3236.    }
  3237.    return(JulialongSetup());
  3238. }
  3239.  
  3240. TrigPlusSqrfpSetup()
  3241. {
  3242.    curfractalspecific->per_pixel =  juliafp_per_pixel;
  3243.    curfractalspecific->orbitcalc =  TrigPlusSqrfpFractal;
  3244.    if(parm.x == 1.0 && parm.y == 0.0 && parm2.y == 0.0 && debugflag != 90)
  3245.    {
  3246.       if(parm2.x == 1.0)    /* Scott variant */
  3247.      curfractalspecific->orbitcalc =  ScottTrigPlusSqrfpFractal;
  3248.       else if(parm2.x == -1.0)    /* Skinner variant */
  3249.      curfractalspecific->orbitcalc =  SkinnerTrigSubSqrfpFractal;
  3250.    }
  3251.    return(JuliafpSetup());
  3252. }
  3253.  
  3254. TrigPlusTriglongSetup()
  3255. {
  3256.    FnPlusFnSym();
  3257.    if(trigndx[1] == SQR)
  3258.       return(TrigPlusSqrlongSetup());
  3259.    curfractalspecific->per_pixel =  long_julia_per_pixel;
  3260.    curfractalspecific->orbitcalc =  TrigPlusTrigFractal;
  3261.    if(lparm.x == fudge && lparm.y == 0L && lparm2.y == 0L && debugflag != 90)
  3262.    {
  3263.       if(lparm2.x == fudge)       /* Scott variant */
  3264.      curfractalspecific->orbitcalc =  ScottTrigPlusTrigFractal;
  3265.       else if(lparm2.x == -fudge)  /* Skinner variant */
  3266.      curfractalspecific->orbitcalc =  SkinnerTrigSubTrigFractal;
  3267.    }
  3268.    return(JulialongSetup());
  3269. }
  3270.  
  3271. TrigPlusTrigfpSetup()
  3272. {
  3273.    FnPlusFnSym();
  3274.    if(trigndx[1] == SQR)
  3275.       return(TrigPlusSqrfpSetup());
  3276.    curfractalspecific->per_pixel =  otherjuliafp_per_pixel;
  3277.    curfractalspecific->orbitcalc =  TrigPlusTrigfpFractal;
  3278.    if(parm.x == 1.0 && parm.y == 0.0 && parm2.y == 0.0 && debugflag != 90)
  3279.    {
  3280.       if(parm2.x == 1.0)    /* Scott variant */
  3281.      curfractalspecific->orbitcalc =  ScottTrigPlusTrigfpFractal;
  3282.       else if(parm2.x == -1.0)    /* Skinner variant */
  3283.      curfractalspecific->orbitcalc =  SkinnerTrigSubTrigfpFractal;
  3284.    }
  3285.    return(JuliafpSetup());
  3286. }
  3287.  
  3288. FnPlusFnSym() /* set symmetry matrix for fn+fn type */
  3289. {
  3290.    static char far fnplusfn[7][7] =
  3291.    {/* fn2 ->sin     cos    sinh    cosh   exp    log    sqr  */
  3292.    /* fn1 */
  3293.    /* sin */ {PI_SYM,XAXIS, XYAXIS, XAXIS, XAXIS, XAXIS, XAXIS},
  3294.    /* cos */ {XAXIS, PI_SYM,XAXIS,  XYAXIS,XAXIS, XAXIS, XAXIS},
  3295.    /* sinh*/ {XYAXIS,XAXIS, XYAXIS, XAXIS, XAXIS, XAXIS, XAXIS},
  3296.    /* cosh*/ {XAXIS, XYAXIS,XAXIS,  XYAXIS,XAXIS, XAXIS, XAXIS},
  3297.    /* exp */ {XAXIS, XYAXIS,XAXIS,  XAXIS, XYAXIS,XAXIS, XAXIS},
  3298.    /* log */ {XAXIS, XAXIS, XAXIS,  XAXIS, XAXIS, XAXIS, XAXIS},
  3299.    /* sqr */ {XAXIS, XAXIS, XAXIS,  XAXIS, XAXIS, XAXIS, XYAXIS}
  3300.    };
  3301.    if(parm.y == 0.0 && parm2.y == 0.0)
  3302.     { if(trigndx[0] < 7 && trigndx[1] < 7)  /* bounds of array JCO 5/6/92*/
  3303.         symmetry = fnplusfn[trigndx[0]][trigndx[1]];  /* JCO 5/6/92 */
  3304.     }                 /* defaults to XAXIS symmetry JCO 5/6/92 */
  3305.    else
  3306.       symmetry = NOSYM;
  3307.    return(0);
  3308. }
  3309.  
  3310. LambdaTrigOrTrigSetup()
  3311. {
  3312. /* default symmetry is ORIGIN  JCO 2/29/92 (changed from PI_SYM) */
  3313.    longparm = &lparm; /* added to consolidate code 10/1/92 JCO */
  3314.    floatparm = &parm;
  3315.    if ((trigndx[0] == EXP) || (trigndx[1] == EXP))
  3316.       symmetry = NOSYM; /* JCO 1/9/93 */
  3317.    if ((trigndx[0] == LOG) || (trigndx[1] == LOG))
  3318.       symmetry = XAXIS;
  3319.    get_julia_attractor (0.0, 0.0);    /* an attractor? */
  3320.    return(1);
  3321. }
  3322.  
  3323. JuliaTrigOrTrigSetup()
  3324. {
  3325. /* default symmetry is XAXIS */
  3326.    longparm = &lparm; /* added to consolidate code 10/1/92 JCO */
  3327.    floatparm = &parm;
  3328.    if(parm.y != 0.0)
  3329.      symmetry = NOSYM;
  3330.    get_julia_attractor (0.0, 0.0);    /* an attractor? */
  3331.    return(1);
  3332. }
  3333.  
  3334. ManlamTrigOrTrigSetup()
  3335. { /* psuedo */
  3336. /* default symmetry is XAXIS */
  3337.    longparm = &linit; /* added to consolidate code 10/1/92 JCO */
  3338.    floatparm = &init;
  3339.    if (trigndx[0] == SQR)
  3340.       symmetry = NOSYM;
  3341.    if ((trigndx[0] == LOG) || (trigndx[1] == LOG))
  3342.       symmetry = NOSYM;
  3343.    return(1);
  3344. }
  3345.  
  3346. MandelTrigOrTrigSetup()
  3347. {
  3348. /* default symmetry is XAXIS_NOPARM */
  3349.    longparm = &linit; /* added to consolidate code 10/1/92 JCO */
  3350.    floatparm = &init;
  3351.    if ((trigndx[0] == 14) || (trigndx[1] == 14)) /* FLIP  JCO 5/28/92 */
  3352.       symmetry = NOSYM;
  3353.    return(1);
  3354. }
  3355.  
  3356.  
  3357. ZXTrigPlusZSetup()
  3358. {
  3359. /*   static char far ZXTrigPlusZSym1[] = */
  3360.    /* fn1 ->  sin   cos    sinh  cosh exp   log   sqr */
  3361. /*         {XAXIS,XYAXIS,XAXIS,XYAXIS,XAXIS,NOSYM,XYAXIS}; */
  3362. /*   static char far ZXTrigPlusZSym2[] = */
  3363.    /* fn1 ->  sin   cos    sinh  cosh exp   log   sqr */
  3364. /*         {NOSYM,ORIGIN,NOSYM,ORIGIN,NOSYM,NOSYM,ORIGIN}; */
  3365.  
  3366.    if(param[1] == 0.0 && param[3] == 0.0)
  3367. /*      symmetry = ZXTrigPlusZSym1[trigndx[0]]; */
  3368.    switch(trigndx[0])
  3369.    {
  3370.       case COS:   /* changed to two case statments and made any added */
  3371.       case COSH:  /* functions default to XAXIS symmetry. JCO 5/25/92 */
  3372.       case SQR:
  3373.       case 9:   /* 'real' cos */
  3374.          symmetry = XYAXIS;
  3375.          break;
  3376.       case 14:   /* FLIP  JCO 2/29/92 */
  3377.          symmetry = YAXIS;
  3378.          break;
  3379.       case LOG:
  3380.          symmetry = NOSYM;
  3381.          break;
  3382.       default:
  3383.          symmetry = XAXIS;
  3384.          break;
  3385.       }
  3386.    else
  3387. /*      symmetry = ZXTrigPlusZSym2[trigndx[0]]; */
  3388.    switch(trigndx[0])
  3389.    {
  3390.       case COS:
  3391.       case COSH:
  3392.       case SQR:
  3393.       case 9:   /* 'real' cos */
  3394.          symmetry = ORIGIN;
  3395.          break;
  3396.       case 14:   /* FLIP  JCO 2/29/92 */
  3397.          symmetry = NOSYM;
  3398.          break;
  3399.       default:
  3400.          symmetry = XAXIS;
  3401.          break;
  3402.       }
  3403.    if(curfractalspecific->isinteger)
  3404.    {
  3405.       curfractalspecific->orbitcalc =  ZXTrigPlusZFractal;
  3406.       if(lparm.x == fudge && lparm.y == 0L && lparm2.y == 0L && debugflag != 90)
  3407.       {
  3408.      if(lparm2.x == fudge)       /* Scott variant */
  3409.          curfractalspecific->orbitcalc =  ScottZXTrigPlusZFractal;
  3410.      else if(lparm2.x == -fudge)  /* Skinner variant */
  3411.          curfractalspecific->orbitcalc =  SkinnerZXTrigSubZFractal;
  3412.       }
  3413.       return(JulialongSetup());
  3414.    }
  3415.    else
  3416.    {
  3417.       curfractalspecific->orbitcalc =  ZXTrigPlusZfpFractal;
  3418.       if(parm.x == 1.0 && parm.y == 0.0 && parm2.y == 0.0 && debugflag != 90)
  3419.       {
  3420.      if(parm2.x == 1.0)    /* Scott variant */
  3421.          curfractalspecific->orbitcalc =  ScottZXTrigPlusZfpFractal;
  3422.      else if(parm2.x == -1.0)    /* Skinner variant */
  3423.          curfractalspecific->orbitcalc =  SkinnerZXTrigSubZfpFractal;
  3424.       }
  3425.    }
  3426.    return(JuliafpSetup());
  3427. }
  3428.  
  3429. LambdaTrigSetup()
  3430. {
  3431.    int isinteger;
  3432.    if((isinteger = curfractalspecific->isinteger))
  3433.       curfractalspecific->orbitcalc =  LambdaTrigFractal;
  3434.    else
  3435.       curfractalspecific->orbitcalc =  LambdaTrigfpFractal;
  3436.    switch(trigndx[0])
  3437.    {
  3438.    case SIN:
  3439.    case COS:
  3440.    case 9:   /* 'real' cos, added this and default for additional functions */
  3441.       symmetry = PI_SYM;
  3442.       if(isinteger)
  3443.      curfractalspecific->orbitcalc =  LambdaTrigFractal1;
  3444.       else
  3445.      curfractalspecific->orbitcalc =  LambdaTrigfpFractal1;
  3446.       break;
  3447.    case SINH:
  3448.    case COSH:
  3449.       symmetry = ORIGIN;
  3450.       if(isinteger)
  3451.      curfractalspecific->orbitcalc =  LambdaTrigFractal2;
  3452.       else
  3453.      curfractalspecific->orbitcalc =  LambdaTrigfpFractal2;
  3454.       break;
  3455.    case SQR:
  3456.       symmetry = ORIGIN;
  3457.       break;
  3458.    case EXP:
  3459.       if(isinteger)
  3460.      curfractalspecific->orbitcalc =  LongLambdaexponentFractal;
  3461.       else
  3462.      curfractalspecific->orbitcalc =  LambdaexponentFractal;
  3463.       symmetry = NOSYM; /* JCO 1/9/93 */
  3464.       break;
  3465.    case LOG:
  3466.       symmetry = NOSYM;
  3467.       break;
  3468.    default:   /* default for additional functions */
  3469.       symmetry = ORIGIN;  /* JCO 5/8/92 */
  3470.       break;
  3471.    }
  3472.    get_julia_attractor (0.0, 0.0);    /* an attractor? */
  3473.    if(isinteger)
  3474.       return(JulialongSetup());
  3475.    else
  3476.       return(JuliafpSetup());
  3477. }
  3478.  
  3479. JuliafnPlusZsqrdSetup()
  3480. {
  3481. /*   static char far fnpluszsqrd[] = */
  3482.    /* fn1 ->  sin   cos    sinh  cosh    sqr    exp   log  */
  3483.    /* sin    {NOSYM,ORIGIN,NOSYM,ORIGIN,ORIGIN,NOSYM,NOSYM}; */
  3484.  
  3485. /*   symmetry = fnpluszsqrd[trigndx[0]];   JCO  5/8/92 */
  3486.    switch(trigndx[0]) /* fix sqr symmetry & add additional functions */
  3487.    {
  3488.    case COS: /* cosxx */
  3489.    case COSH:
  3490.    case SQR:
  3491.    case 9:   /* 'real' cos */
  3492.    case 10:  /* tan */
  3493.    case 11:  /* tanh */
  3494.    symmetry = ORIGIN;
  3495.     /* default is for NOSYM symmetry */
  3496.    }
  3497.    if(curfractalspecific->isinteger)
  3498.       return(JulialongSetup());
  3499.    else
  3500.       return(JuliafpSetup());
  3501. }
  3502.  
  3503. SqrTrigSetup()
  3504. {
  3505. /*   static char far SqrTrigSym[] = */
  3506.    /* fn1 ->  sin    cos    sinh   cosh   sqr     exp   log  */
  3507. /*         {PI_SYM,PI_SYM,XYAXIS,XYAXIS,XYAXIS,XAXIS,XAXIS}; */
  3508. /*   symmetry = SqrTrigSym[trigndx[0]];      JCO  5/9/92 */
  3509.    switch(trigndx[0]) /* fix sqr symmetry & add additional functions */
  3510.    {
  3511.    case SIN:
  3512.    case COS: /* cosxx */
  3513.    case 9:   /* 'real' cos */
  3514.    symmetry = PI_SYM;
  3515.     /* default is for XAXIS symmetry */
  3516.    }
  3517.    if(curfractalspecific->isinteger)
  3518.       return(JulialongSetup());
  3519.    else
  3520.       return(JuliafpSetup());
  3521. }
  3522.  
  3523. FnXFnSetup()
  3524. {
  3525.    static char far fnxfn[7][7] =
  3526.    {/* fn2 ->sin     cos    sinh    cosh  exp    log    sqr */
  3527.    /* fn1 */
  3528.    /* sin */ {PI_SYM,YAXIS, XYAXIS,XYAXIS,XAXIS, NOSYM, XYAXIS},
  3529.    /* cos */ {YAXIS, PI_SYM,XYAXIS,XYAXIS,XAXIS, NOSYM, XYAXIS},
  3530.    /* sinh*/ {XYAXIS,XYAXIS,XYAXIS,XYAXIS,XAXIS, NOSYM, XYAXIS},
  3531.    /* cosh*/ {XYAXIS,XYAXIS,XYAXIS,XYAXIS,XAXIS, NOSYM, XYAXIS},
  3532.    /* exp */ {XAXIS, XAXIS, XAXIS, XAXIS, XAXIS, NOSYM, XYAXIS},
  3533.    /* log */ {NOSYM, NOSYM, NOSYM, NOSYM, NOSYM, XAXIS, NOSYM},
  3534.    /* sqr */ {XYAXIS,XYAXIS,XYAXIS,XYAXIS,XYAXIS,NOSYM, XYAXIS},
  3535.    };
  3536.    /*
  3537.    if(trigndx[0]==EXP || trigndx[0]==LOG || trigndx[1]==EXP || trigndx[1]==LOG)
  3538.       symmetry = XAXIS;
  3539.    else if((trigndx[0]==SIN && trigndx[1]==SIN)||(trigndx[0]==COS && trigndx[1]==COS))
  3540.       symmetry = PI_SYM;
  3541.    else if((trigndx[0]==SIN && trigndx[1]==COS)||(trigndx[0]==COS && trigndx[1]==SIN))
  3542.       symmetry = YAXIS;
  3543.    else
  3544.       symmetry = XYAXIS;
  3545.    */
  3546.    if(trigndx[0] < 7 && trigndx[1] < 7)  /* bounds of array JCO 5/22/92*/
  3547.         symmetry = fnxfn[trigndx[0]][trigndx[1]];  /* JCO 5/22/92 */
  3548.                     /* defaults to XAXIS symmetry JCO 5/22/92 */
  3549.    else {  /* added to complete the symmetry JCO 5/22/92 */
  3550.       if (trigndx[0]==LOG || trigndx[1] ==LOG) symmetry = NOSYM;
  3551.       if (trigndx[0]==9 || trigndx[1] ==9) { /* 'real' cos */
  3552.          if (trigndx[0]==SIN || trigndx[1] ==SIN) symmetry = PI_SYM;
  3553.          if (trigndx[0]==COS || trigndx[1] ==COS) symmetry = PI_SYM;
  3554.       }
  3555.       if (trigndx[0]==9 && trigndx[1] ==9) symmetry = PI_SYM;
  3556.    }
  3557.    if(curfractalspecific->isinteger)
  3558.       return(JulialongSetup());
  3559.    else
  3560.       return(JuliafpSetup());
  3561. }
  3562.  
  3563. MandelTrigSetup()
  3564. {
  3565.    int isinteger;
  3566.    if((isinteger = curfractalspecific->isinteger))
  3567.       curfractalspecific->orbitcalc =  LambdaTrigFractal;
  3568.    else
  3569.       curfractalspecific->orbitcalc =  LambdaTrigfpFractal;
  3570.    symmetry = XYAXIS_NOPARM;
  3571.    switch(trigndx[0])
  3572.    {
  3573.    case SIN:
  3574.    case COS:
  3575.       if(isinteger)
  3576.      curfractalspecific->orbitcalc =  LambdaTrigFractal1;
  3577.       else
  3578.      curfractalspecific->orbitcalc =  LambdaTrigfpFractal1;
  3579.       break;
  3580.    case SINH:
  3581.    case COSH:
  3582.       if(isinteger)
  3583.      curfractalspecific->orbitcalc =  LambdaTrigFractal2;
  3584.       else
  3585.      curfractalspecific->orbitcalc =  LambdaTrigfpFractal2;
  3586.       break;
  3587.    case EXP:
  3588.       symmetry = XAXIS_NOPARM;
  3589.       if(isinteger)
  3590.      curfractalspecific->orbitcalc =  LongLambdaexponentFractal;
  3591.       else
  3592.      curfractalspecific->orbitcalc =  LambdaexponentFractal;
  3593.       break;
  3594.    case LOG:
  3595.       symmetry = XAXIS_NOPARM;
  3596.       break;
  3597.    default:   /* added for additional functions, JCO 5/25/92 */
  3598.       symmetry = XYAXIS_NOPARM;
  3599.       break;
  3600.    }
  3601.    if(isinteger)
  3602.       return(MandellongSetup());
  3603.    else
  3604.       return(MandelfpSetup());
  3605. }
  3606.  
  3607. MarksJuliaSetup()
  3608. {
  3609.    c_exp = param[2];
  3610.    longparm = &lparm;
  3611.    lold = *longparm;
  3612.    if(c_exp > 2)
  3613.       lcpower(&lold,c_exp,&lcoefficient,bitshift);
  3614.    else if(c_exp == 2)
  3615.    {
  3616.       lcoefficient.x = multiply(lold.x,lold.x,bitshift) - multiply(lold.y,lold.y,bitshift);
  3617.       lcoefficient.y = multiply(lold.x,lold.y,bitshiftless1);
  3618.    }
  3619.    else if(c_exp < 2)
  3620.       lcoefficient = lold;
  3621.    get_julia_attractor (0.0, 0.0);    /* an attractor? */
  3622.    return(1);
  3623. }
  3624.  
  3625. MarksJuliafpSetup()
  3626. {
  3627.    c_exp = param[2];
  3628.    floatparm = &parm;
  3629.    old = *floatparm;
  3630.    if(c_exp > 2)
  3631.       cpower(&old,c_exp,&coefficient);
  3632.    else if(c_exp == 2)
  3633.    {
  3634.       coefficient.x = sqr(old.x) - sqr(old.y);
  3635.       coefficient.y = old.x * old.y * 2;
  3636.    }
  3637.    else if(c_exp < 2)
  3638.       coefficient = old;
  3639.    get_julia_attractor (0.0, 0.0);    /* an attractor? */
  3640.    return(1);
  3641. }
  3642.  
  3643. SierpinskiSetup()
  3644. {
  3645.    /* sierpinski */
  3646.    periodicitycheck = 0;        /* disable periodicity checks */
  3647.    ltmp.x = 1;
  3648.    ltmp.x = ltmp.x << bitshift; /* ltmp.x = 1 */
  3649.    ltmp.y = ltmp.x >> 1;            /* ltmp.y = .5 */
  3650.    return(1);
  3651. }
  3652.  
  3653. SierpinskiFPSetup()
  3654. {
  3655.    /* sierpinski */
  3656.    periodicitycheck = 0;        /* disable periodicity checks */
  3657.    tmp.x = 1;
  3658.    tmp.y = 0.5;
  3659.    return(1);
  3660. }
  3661.  
  3662. extern char usr_floatflag;
  3663.  
  3664. HalleySetup()
  3665. {
  3666.    /* Halley */
  3667.    periodicitycheck=0;
  3668.  
  3669.    if(usr_floatflag)
  3670.      fractype = HALLEY; /* float on */
  3671.    else
  3672.      fractype = MPHALLEY;
  3673.  
  3674.    curfractalspecific = &fractalspecific[fractype];
  3675.  
  3676.    degree = (int)parm.x;
  3677.    if(degree < 2)
  3678.       degree = 2;
  3679.    param[0] = (double)degree;
  3680.  
  3681. /*  precalculated values */
  3682.    AplusOne = degree + 1; /* a+1 */
  3683.    Ap1deg = AplusOne * degree;
  3684.  
  3685. #ifndef XFRACT
  3686.    if(fractype == MPHALLEY) {
  3687.       setMPfunctions();
  3688.       mpAplusOne = *pd2MP((double)AplusOne);
  3689.       mpAp1deg = *pd2MP((double)Ap1deg);
  3690.       mptmpparmy = *pd2MP(parm.y);
  3691.       mptmpparm2x = *pd2MP(parm2.x);
  3692.       mpone       = *pd2MP(1.0);
  3693.    }
  3694. #endif
  3695.  
  3696.    if(degree % 2)
  3697.      symmetry = XAXIS;   /* odd */
  3698.    else
  3699.      symmetry = XYAXIS; /* even */
  3700.    return(1);
  3701. }
  3702.  
  3703. PhoenixSetup()
  3704. {
  3705.    longparm = &lparm; /* added to consolidate code 10/1/92 JCO */
  3706.    floatparm = &parm;
  3707.    degree = (int)parm2.x;
  3708.    if(degree < 2 && degree > -3) degree = 0;
  3709.    param[2] = (double)degree;
  3710.    if(degree == 0){
  3711.      if(usr_floatflag)
  3712.        curfractalspecific->orbitcalc =  PhoenixFractal;
  3713.      else
  3714.        curfractalspecific->orbitcalc =  LongPhoenixFractal;
  3715.    }
  3716.    if(degree >= 2){
  3717.      degree = degree - 1;
  3718.      if(usr_floatflag)
  3719.        curfractalspecific->orbitcalc =  PhoenixPlusFractal;
  3720.      else
  3721.        curfractalspecific->orbitcalc =  LongPhoenixPlusFractal;
  3722.    }
  3723.    if(degree <= -3){
  3724.      degree = abs(degree) - 2;
  3725.      if(degree % 2)
  3726.        symmetry = XYAXIS; /* odd */
  3727.      else
  3728.        symmetry = XAXIS; /* even */
  3729.      if(usr_floatflag)
  3730.        curfractalspecific->orbitcalc =  PhoenixMinusFractal;
  3731.      else
  3732.        curfractalspecific->orbitcalc =  LongPhoenixMinusFractal;
  3733.    }
  3734.  
  3735.    return(1);
  3736. }
  3737.  
  3738. MandPhoenixSetup()
  3739. {
  3740.    longparm = &linit; /* added to consolidate code 10/1/92 JCO */
  3741.    floatparm = &init;
  3742.    degree = (int)parm2.x;
  3743.    if(degree < 2 && degree > -3) degree = 0;
  3744.    param[2] = (double)degree;
  3745.    if(degree == 0){
  3746.      if(usr_floatflag)
  3747.        curfractalspecific->orbitcalc =  PhoenixFractal;
  3748.      else
  3749.        curfractalspecific->orbitcalc =  LongPhoenixFractal;
  3750.    }
  3751.    if(degree >= 2){
  3752.      degree = degree - 1;
  3753.      if(usr_floatflag)
  3754.        curfractalspecific->orbitcalc =  PhoenixPlusFractal;
  3755.      else
  3756.        curfractalspecific->orbitcalc =  LongPhoenixPlusFractal;
  3757.    }
  3758.    if(degree <= -3){
  3759.      degree = abs(degree) - 2;
  3760.      if(usr_floatflag)
  3761.        curfractalspecific->orbitcalc =  PhoenixMinusFractal;
  3762.      else
  3763.        curfractalspecific->orbitcalc =  LongPhoenixMinusFractal;
  3764.    }
  3765.  
  3766.    return(1);
  3767. }
  3768.  
  3769. StandardSetup()
  3770. {
  3771.    if(fractype==UNITYFP)
  3772.       periodicitycheck=0;
  3773.    return(1);
  3774. }
  3775.  
  3776. demowalk()
  3777. {
  3778.     extern double param[];        /* optional user parameters */
  3779.     extern int maxit;            /* maximum iterations (steps) */
  3780.     extern int rflag, rseed;        /* random number seed */
  3781.     extern int xdots, ydots;        /* image coordinates */
  3782.     extern int colors;                  /* maximum colors available */
  3783.     extern double far *dx0, far *dy0;    /* arrays of pixel coordinates */
  3784.     extern double far *dx1, far *dy1;    /* (... for skewed zoom-boxes) */
  3785.  
  3786.     float stepsize;            /* average stepsize */
  3787.     int xwalk, ywalk;            /* current position */
  3788.     int xstep, ystep;            /* current step */
  3789.     int steps;                /* number of steps */
  3790.     int color;                /* color to draw this step */
  3791.     float temp, tempadjust;        /* temporary variables */
  3792.  
  3793. if (param[0] != 999) {            /* if 999, do a Mandelbrot instead */
  3794.  
  3795.     srand(rseed);            /* seed the random number generator */
  3796.     if (!rflag) ++rseed;
  3797.     tempadjust = RAND_MAX >> 2;        /* adjustment factor */
  3798.  
  3799.     xwalk = xdots / 2;            /* start in the center of the image */
  3800.     ywalk = ydots / 2;
  3801.  
  3802.     stepsize = min(xdots, ydots)     /* calculate average stepsize */
  3803.                * (param[0]/100.0);    /* as a percentage of the image */
  3804.  
  3805.     color = max(0, min(colors, param[1]));  /* set the initial color */  
  3806.  
  3807.     for (steps = 0; steps < maxit; steps++) { /* take maxit steps */
  3808.         if (keypressed())        /* abort if told to do so */
  3809.             return(0);
  3810.         temp = rand();            /* calculate the next xstep */
  3811.         xstep = ((temp/tempadjust) - 2.0) * stepsize;
  3812.         xstep = min(xwalk + xstep, xdots - 1);
  3813.         xstep = max(0, xstep);
  3814.         temp = rand();            /* calculate the next ystep */
  3815.         ystep = ((temp/tempadjust) - 2.0) * stepsize;
  3816.         ystep = min(ywalk + ystep, ydots - 1);
  3817.         ystep = max(0, ystep);
  3818.         if (param[1] == 0.0)        /* rotate the colors? */
  3819.             if (++color >= colors)    /* rotate the colors, avoiding */
  3820.                 color = 1;        /* the background color 0 */
  3821.         /* the draw_line function is borrowed from the 3D routines */
  3822.         draw_line(xwalk, ywalk,xstep,ystep,color);
  3823.         /* or, we could be on a pogo stick and just displaying
  3824.            where we landed...
  3825.         putcolor(xstep, ystep, color);
  3826.         */
  3827.  
  3828.         xwalk = xstep;            /* remember where we were */
  3829.         ywalk = ystep;
  3830.         }
  3831.     return(1);                          /* we're done */
  3832.  
  3833. } else {                /* a simple Mandelbrot routine */
  3834.  
  3835.     /* the following routine determines the X and Y values of
  3836.        each pixel coordinate and calculates a simple mandelbrot
  3837.        fractal with them - slowly, but surely */
  3838.     int ix, iy;
  3839.     for (iy = 0; iy < ydots; iy++) {
  3840.         for (ix = 0; ix < xdots; ix++) {
  3841.             int iter;
  3842.             double x, y, newx, newy, tempxx, tempxy, tempyy;
  3843.             /* first, obtain the X and Y coordinate values of this pixel */
  3844.             x = dx0[ix]+dx1[iy];
  3845.             y = dy0[iy]+dy1[ix];
  3846.             /* now initialize the temporary values */
  3847.             tempxx = tempyy = tempxy = 0.0;
  3848.             if (keypressed())        /* abort if told to do so */
  3849.                 return(0);
  3850.             /* the inner iteration loop */
  3851.             for (iter = 1; iter < maxit; iter++) {
  3852.                 /* calculate the X and Y values of Z(iter) */
  3853.                 newx = tempxx - tempyy + x;
  3854.                 newy = tempxy + tempxy + y;
  3855.                 /* calculate the temporary values */
  3856.                 tempxx = newx * newx;
  3857.                 tempyy = newy * newy;
  3858.                 tempxy = newx * newy;
  3859.                 /* are we done yet? */
  3860.                 if (tempxx + tempyy > 4.0) break;
  3861.                 }
  3862.             /* color in the pixel */
  3863.             putcolor(ix, iy, iter & (colors - 1));
  3864.             }
  3865.         }
  3866.     return(1);                          /* we're done */
  3867.     }
  3868.  
  3869. }
  3870.